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ABSTRACT 


This thesis investigates continucus-spatial dynamic models of 
the composition behavior of binary plate distillation columns, (‘These 
partial differential equation models are developed from basic distil- 
lation column principles and the discreteeplate models by treating 
the plate number as a continuousespatial variable. Linearized 
continuous=spatial models are investigated in detail. 


The central purpose of this thesis is the development, pres= 
entation, suggested analytical solution technique, and example column 
evaluation of the Linear Polynomial-Coefficient Model (LPCM) which 
is a linearized continuousespatial model in which the coefficients 
in the partial differential equation are n-th degree polynomials in 
the spatial variable, A general analytical solution technique is 
proposed in which the spatial differential eigenvalue problem 
resulting from separation of variables is transformed to a Liouville 
Normal-Form equation which is then converted to a homogeneous 
Fredholm II integral equation, Several simple examples of the model 
are soived in complete detail, and the proposed solution technique 
is applied to a model with first-degree polynomial coefficients. 

An analytical and a computational analysis giving the details of 
each step of the solution technique is presented, It is suggested 
that greatly reduced computation times compared to discrete models 
will result from application of the proposed solution technique. 


A bibliography of 352 references, 202 of wnich pertain directly 
to distillation column dynamics and control, is presented and related 
to the areas of the thesis, 


Thesis Supervisor: Lawrence B, Evans 
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SECTION 1 
INTRODUCTION TO DISTILLATION (1) 
Tl BASIC PRINCIPL2S OF BINARY DISTILLATION 
I2 THE BINARY PLATE DISTILLATION COLUMN 
I3 PHILOSOPHY AND DISTILLATION 


T4 A BRIEF REVIEW OF THE LITERATURE OF DISTILLATION 


THE FUNDAMENTAL AND UNIVERSAL ABSOLUTE: 

"EXISTENCE, REALITY, THE EXTERNAL WORLD, IS WHAT IT IS, INDEPEN- 
DENT OF MAN*S CONSCIOUSNESS, INDEPENDENT OF ANYONE'S KNOWLEDGE, 
JUDGMENT, BELIEFS, HOPES, WISHES, OR FEARS — THAT FACTS ARE FACTS, 
THAT A IS A, THAT THINGS ARE WHAT THEY ARE." 


NATHANIEL BRANDEN (B-28) 


THIS SECTION PRESENTS THE BASIC PHYSICAL AND MATHEMATICAL IDEAS 
OF BINARY DISTILLATION ORIENTED TOWARD FORMING A MATHEMATICAL MODEL 
OF A BINARY PLATE DISTILLATION COLUMN. ‘THE GENERAL PROBLEM OF OBSERV- 
ING AND MODELING PHYSICAL SYSTEMS IS DISCUSSED AND THE HISTORY AND 
GENERAL LITERATURE OF DISTILLATION ARE REVIEWED BRIEFLY. SEVERAL 
OPINIONS ON THE PHILOSOPHY OF OBSERVING REALITY AS RELATED TO MODELING 
A DISTILLATION COLUMN AND ON THE DESIRE TO ACHIEVE PROFIT AS RELATED 
TO THE COST OF SEPARATION OF COMPONENTS ARE PRESENTED, 

THE CENTRAL PURPOSE OF THIS THESIS IS THE DEVELOPMENT, PRESENTATION, 
SUGGESTED SOLUTION TECHNIQUE, AND EVALUATION OF THE LINEAR POLYNOMIAL- 
COEFFICIENT MODEL (LPCM) OF THE DYNAMIC BEHAVIOR OF A BINARY PLATE 
DISTILLATION COLUMN, THIS SECTION PRESENTS THE BASIC PRINCIPLES LEADING 
UP TO THE COMPLETE DEVELOPMENT OF THE LPCM IN SECTION 2(M) AND THE 


PRESENTATION OF AN INTEGRAL EQUATION SOLUTION TECHNIQUE IN SECTION3(L). 
ale 





CHAPTER I1 


BASIC PRINCIPLES OF BINARY DISTILLATION 

The purpose of Chapters Il and I2 is to present the basic 
principles behind the separation of binary mixtures and to develop 
the techniques for describing mathematically and graphically the 
steady state characteristics of a distillation column. These topics 
are presented for several reasons: 

le Some readers may not be acquainted with the basic theory, 
which is necessary for the developments which follow. 

2e This is a convenient method for developing a consistent 
notation for use throughout this thesis. 

3. Presentations of those basic principles specifically 
used in this thesis are not often found in the existing literature 
in forms easily understood by readers without some background in 
the subject. 

There are many available techniques for describing binary 
distillation in quantitative or qualitative terms. Most of these 
techniques are oriented specifically toward one of these two outlooks. 
The technique to be described in Chapters Il and I2 is the use of the 
McCabe-Thiele diagram which has the unique advantage of a quantitative, 
qualitative, and, in a sense, visual insight into binary distillation. 

The theory developed in this chapter represents a combination 
of the developments available in the literature. The textbooks used 
in this presentation of the theory in this chapter are listed below 
in the order of decreasing pertinence to this chapter, 


Gould (G~3) 


Sie 





Bennett and Myers (B-1) 
Van Winkle (V-1) 
Holland (H-2), (H-7) 
Some of the other textbooks which present the basic principles of 
distillation are: (T-1), (S-1), (H=-1), (C-1), (R-17), (c-5), (H-6), 
(H=5), (R-21), and (M-2). 
Il.l1 DEFINITION OF BINARY DISTILLATION AND RELATIVE VOLATILITY 
Binary distillation is defined as a process which separates a 
mixture of two components by utilizing mass transfer between the 
liquid and vapor phases of the components. The essence of this sep- 
aration lies in the fact that when the vapor and liquid phases of a 
binary mixture are in equilibrium, the vapor is richer in the lighter 
component than is the liquid. The process by which the vapor phase 
becomes richer in the component which boils at the lower temperature 
(lighter component) is called mass transfer. Although the term 
distillation is occasionally used to describe the removal of volatile 
materials from solids, the term as used in this thesis will apply only 
to the separation of volatile components found in liquid solutions. 
The equilibrium mentioned above is a phase equilibrium in which 
the properties of the two phases depend upon the physical characteristics 
of the two components which are present. The primary physical char- 
acteristic of interest in this thesis is the relative volatility q. 
If the concentration of the lighter component in the liquid phase is 
defined as u and the concentration of the lighter component in the 
vapor phase is defined as f(u), then the relative volatility is given 


by equation I1.1. This expression is only valid when q does not vary 


eS 





azwi-u, flu 1 ee 
u le-f(u 


with composition, which is a valid approximation for a large number 
of liquids. Equation Il].1, then, defines equilibrium between the two 
phases in terms of the constant relative volatility aq. 
I1.2 THE EQUILIBRIUM CURVE 

The assumption that q is constant is essential to the develop- 
ments of this thesis because it allows a specific equilibrium function 
to represent the characteristics of the phase mixture. Thus, an 
equilibrium phase diagram can be drawn as in Figure I1.1 and a specific 
equilibrium function f(u) can be expressed in Equation 11.2, 


f(u) = au i ye 
1+ (qe )u 
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Figure T1.1 - THE EQUILIBRIUM CURVE 








In Figure Il.1 the curved line (equilibrium curve) represents 
a graphical plot of equation 11.2. If the relative volatility is not 
assumed to be constant, then this curve may have a significantly 
different shape and be represented by a different functional relation- 
ship. In trying to understand the meaning of the equilibrium curve 
it is helpful to consider a point u, in Figure I1.1. 

Given a binary mixture which has liquid and vapor phases at 
equilibrium, then the concentration of the lighter component in the 
liquid is given by Us» and the concentration of the lighter component 
in the vapor is given by f(u,). Thus, going from the liquid to the 
vapor in the two phase mixture corresponds to going from the u line 
to the f(u) line on the equilibrium curve. The fact that the relative 
volatility is such that f(u) > u in this case implies that a partial 
separation of the mixture can be accomplished by separating the vapor 
from the liquid. 

T1.3 USING THE EQUILIBRIUM CURVE TO DESCRIBE SEPARATION 

If the vapor is separated from the liquid in a two phase binary 
mixture, then condensing the vapor produces two separate liquids of 
different compositions, the one having been condensed being richer in 
the lighter component. This method of separation is the key to binary 
distillation. 

This separation method can be visualized by using the equilibrium 
curve in Figure I1.2. The original liquid at composition u, is boiled 
and part of it becomes vapor at composition f(u,). If the vapor at 
f(u,) is then separated from the liquid at u, and then condensed, the 


liquid condensate is of composition u, = f(u,). The new liquid at 


eye 





composition u, can now be boiled to produce vapor of composition 


f(u,), where £(u, ) >u,. Thus, by successively boiling and condensing 
the liquids and vapors, a desired degree of separation in terms of 


the upper and lower (on the graph of Figure I1.2) liquids can be achieved. 
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Figure 11.2 + DESCRIBING SEPARATION ON THE EQUILIBRIUM CURVE 
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T1.4 PHYSICAL SEPARATION OF VAPOR AND LIQUID 
The next problem to be considered is that of the physical equip- 


ment necessary to separate the vapor of f(u,) from the liquid at u., 
while at the es time allowing the liquid-vapor contact between the 
vapor at f(u,) and the liquid at u,- There are two commonly used 
types of equipment for accomplishing this: packed columns and plate 
columns, both of which have innumerable variations. 

Packed columns use a form of continuous contacting of liquid 
and vapor phases. This is accomplished by use of many small devices 
in the form of rings, saddles, spheres, etc,, which are randomly 
"packed" into the column. This arrangement is designed to provide a 
very large surface area for a given tower volume because the contacting 
and mass transfer occur at the surfaces of the particles, The 
mathematical models (Section M) of packed columns are usually partial- 
differential equations. Packed columns will not be investigated 
specifically in this thesis, but it is expected that the developments 
of this thesis could be applied to packed columns, 

Plate columns accomplish vapor-liquid separation and vapor-liquid 
contacting by allowing the liquid at composition uw to flow over the 
top of a plate, which is designed to let the vapor at f(u,) pass 
through into the liquid flowing over it, There are a multitude of 
different plate designs which accomplish this efficiently, such as 
bubble-cap, perforated, and sieve plates. Figure I1.3 shows the basic 
physical characteristics of a bubble-cap plate. The concentrations 
shown in Figure I1.3 correspond to those on the equilibrium curve of 


Figure I1.2. Using these two figures one can see both physically and 
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conceptually how separation is achieved. 

Physically, the liquid at u, enters the plate of Figure Il1.3 from 
the left-side downcomer of the plate above, The Uy liquid then flows 
across the plate, in this case left to right, and is mixed (contacted) 
with the vapor f(u,) coming from beneath the plate. These two phases 
reach an equilibrium such that the liquid leaving the plate is at u, 
and the vapor leaving the liquid is at f(u,). Since f(uj) > uy, a 
partial separation of the components results. Conceptually, following 
these same arguments on the lines in Figure I1.2 reveals the properties 


of the mechanism of separation in terms of the mathematical expressions 
=o0= 





describing equilibrium. 

This chapter has presented the basic principles of binary 
separation using distillation. This separation has been described 
conceptually in terms of the relative volatility and the resulting 
equilibrium curve and physically in terms of the operation of a 
bpubble-cap plate. The next chapter puts several of these plates 
together and presents the physical and conceptual descriptions of 


a complete distillation column. 
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CHAPTER I2 


THE BINARY PLATE DISTILLATION COLUMN 

The particular binary separation device chosen for study in this 
thesis is the plate distillation column. A plate distillation column 
is a vertical cascaded arrangement of individual plates (Figure I1.3) 
supported in a tower or column. The distillation column alone cannot 
perform separation but requires auxiliary equipment for its operation. 
This chapter presents a very brief description of the physical operation 
of a distillation column and its auxiliary equipment, a graphical 
description of its operation in terms of the McCabe-Thiele diagram, 
and a steady-state mathematical model of its operation. 

I2.1 PHYSICAL OPERATION OF A BINARY PLATE DISTILLATION COLUMN 

The physical separation of a binary mixture by an individual plate 
Was described in Chapter Il. In general, one plate is usually in- 
adequate to achieve the desired degree of separation of the mixture, 
thus many plates are combined in cascade to achieve greater output 
purity, An arrangement of eleven such plates is shown in the column 
of Figure 12,1. On any given plate in this column the liquid and vapor 
mix and reach equilibrium, with the vapor rising in the column plate- 
by-plate through the bubble caps and with the liquid flowing back and 
forth down the column. 

Two of the auxiliary equipments necessary to the operation of the 
column are the condenser and the reboiler. The vapor flowing out of 
the top of the column enters the condenser where it is liquified and the 
resulting liquid is partially fed back to the top tray and partially 
removed as tops product, Similarly, the liquid flowing out the bottom 
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of the column is divided into that portion removed as bottoms product 
and that portion which is vaporized in the reboiler and fed back to the 
bottom plate, The other auxiliary equipments, such as pumps to cir- 
culate the fluid, valves for control, structural supports, and many 
others, are not shown in Figure [2,1] but are, nevertheless, essential 
to the operation of the systen. 

Having described the operation of the internal cycling of the 
fluids in the column, the overall operation of the system can now be 
discussed. The input mixture to be separated is fed to that tray 
designated as the feed tray and enters the fluid cycle. The separated 
outputs are then taken off as the top and bottom products. Energy for 
operation of the system is supplied by the reboiler, steam heated in 
the case of the column of Figure 12.1, and energy is removed from the 
system by the condenser, water cooled in Figure I2.l. 

There are innumerable types, arrangements, designs, and sizes of 
distillation columns, just as there are many different aspects of any 
given column which can be studied, such as chemical, thermal, structural, 
fluidic, economic, and environmental serie: The central purpose of 
any column is to separate a binary mixture and the main criterion of 
“goodness” of the operation of the column is how well it performs this 
separation. Distillation columns are usually designed to separate the 
mixture to desired purity in such a way as to maximize the economic 
profit derived from the sale of the separated products. 

Output variations in the product purity greatly affect the profits 
derived from such sales, If the output is overly pure, then the revenue 


resulting from the sale of material thought to be of lower purity will 
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not be as high as it could have been, If the output is lower than the 
desired composition, then it cannot be sold at the price set for 
material of higher quality and profits will be reduced, Thus, composi- 
tion variations in distillation column inputs and outputs are rather 
important, and mathematical models for the interaction between feed 
composition changes and output composition changes will be developed 

in this thesis for the purpose of studying such variations. 

12.2 THE MCCABE-THIELE DIAGRAM 

One of the best graphical techniques for describing the steady- 
state operation of a binary plate column is the McCabe-Thiele diagram 
(M-10). A technique for describing the steady-state operation of the 
column is important to considerations involving variations in composi- 
tions, especially when composition variations are to be interpreted 
as upsets or deviations from an initial steady state to a final steady 
state, The basis of the McCabe-Thiele diagram is the equilibrium curve 
presented in Figure I1.2 for the operation of one plate. For the eleven 
(11) plate column of Figure 12,1 the McCabe-Thiele diagram for gq = 3.0 
(M-15) is shown in Figure 12.2. 

The graphical expression of the individual plate compositions 
presented by the McCabe-Thiele diagram gives both a qualitative and a 
quantitative view of the steady-state operation of the column, Qualita- 
tively, it shows how the individual plates act to increase the purities 
of the output streams by operating in cascade with each step up the 
diagram representing a plate corresponding to a step up the actual 
column. Quantitatively, the individual plate steady-state compositions 


and the operating line slopes can be read directly from the diagran. 
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The actual operating characteristics of the column are evident 
from the McCabe-Thiele Diagram. The amounts or percentages of liquid 
and vapor reflux flows determine the slopes of the two operating lines, 
the upper and lower lines in Figure 12.2. Each point of intersection 
on the operating lines represents the liquid composition of a plate, 
and thus the number of liquid points equals the number of plates in 
the column, lkach point of intersection on the equilibrium curve repre- 
sents the vapor composition between two of the column plates, The 
intersection of the two operating lines is determined by the q-line 
which depends upon the properties and condition of the feed at Ups All 
of these properties are represented by quantities and equations in the 
steady-state model of the distiliation column, 

12.3 A STEADY-STATE MODEL OF A BINARY PLATE DISTILLATION COLUMN 

This subsection begins the developments leading to mathematical 
equations describing the operation of a binary plate distillation column. 
A mathematical model of the steady-state operation is to be presented 
and explained. Since most of the detailed developments of steady-state 
models are well covered in the literature(See B2.1 = Steady-State 
Analysis and McCabe-Thiele Diagrams), the discrete-plate steady-state 
model and its symbols will mereiy be presented, not derived in detail, 
in this section. The format of this presentation will be to present 
the description of the symbols to be used, present the model, and give 
a brief explanation of the model and the assumptions behind it, 

The presentation of the discrete-plate steady-state model begins 
by listing and describing the mathematical symbols to be used in the 


model, This list of symbols is found ir Table I2.1. The symbols are 
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f(u) - 


E, = f'(u,) ~ Un ey J 
Cae 


Where f*(u,) is a non-equilibrium state 


Concentration of the lighter component in the liquid (lbm 
mole lighter comp./lbm liquid) 


Concentration of the lighter component in the vapor (1bm 
mole lighter comp./lbm vapor) 


Relative volatility (dimensionless); assumed constant 
Feed rate (lbm liquid/hour) 

Distillate rate (lbm liquid/hour) 

Withdrawal rate of Bottoms Product (lbm liquid/hour) 


Liquid rate in the upper section (lbm liquid/hour); assumed 
constant 


Liquid rate in the lower section (lbm liquid/hour) ; assumed 
constant 


Vapor rate in the column (lbm vapor/hour); assumed constant 
Upper reflux ratio (lbm liquid/lbm vapor), L,/V 

Lower reflux ratio (lbm liquid/lbm vapor), L,/V 

Feed plate index (number); integer 

Internal plate index (number); integer l en <N 


Portion of the feed which adds to the lower liquid rate; 


L, = L + ran 


Table [2.1 + LIST OF SYMBOLS USD IN DISTILLATION 
COLUMN STEADY-STATE MODELS 
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fairly standard throughout the literature with one important exception. 
Usually, in the literature, the composition of the lighter component 

in the liquid is represented by x and that in the vapor by y. In this 
thesis x will be used in Sections M and L to represent a continuous- 
spatial variable or continuous-plate-number variable; and y is not used, 
but the term f(u) is used instead. The term f(u) is standard in most 
of the literature, The symbols in Table I2.1 also apply to the corres- 
ponding symbols placed on Figure I[2.1. 

Figure I2.3 presents the discrete-plate steady-state model of a 
binary plate distillation column. This model is merely a combination 
of the conservation-of-mass equations (Inflow = Outflow) for each plate 
and a listing of the end, feed, and equilibrium description. These 
equations say mathematically what the McCabe-Thiele diagram says graph- 
ically. Solving these equations mathematically for the plate composi- 
tions is also equivalent to “stepping off" the compositions on the 
McCabe-Thiele diagram. Tne amount of mathematical manipulation involved 
in solving for the plate compositions when there are many trays is 
rather large when compared with the ease of stepping off these numbers 
on the McCabe-Thiele diagran. 

Several very significant assumptions are implicit in the model 
statement of Figure 12.3. The upper and lower liquid rates, the 
column. vapor rate, and the relative volatility have all been assumed 
to remain constant. In addition, the Murphree (M-8) efficiencies E. 
defined in equation I2.1 nave all been assumed to be unity, or equiv=- 
alently, the assumption is made that each plate operates at complete 


equilibrium, 
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Model Equations 
Rectification Section (Upper) 


f(u__s) = Bu + (1-3, )u k+leneN 


d 


Stripping Section (Lower) 
f(u_,) = Buu, + (1-3, )u, lenek 


Where B= L/V and B, = L,/V 


Equilibrium 


qu, 
Eee) 


~ 1+(q-1)u, 


End Conditions 


f(uy) = Ug Lop n=WN 
W = Up Feed n=k, q=l 
uw = U. Bottom n=l 


Feed Condition 


Peet ae or 8 = 8 + aF/V ; q = 1 in this case. 


~ DISCRETE+PLATE STEADY-STATE MODEL 





To summarize briefly, Chapters Ii and I2 present the basic 
principles and mathematical descriptions for the operation of a binary 
plate distillation column. This presentation takes the form of a 
graphical description in terms of equilibrium curves and the McCabe- 
Thiele diagram, a brief physical description, and a set of mathematical 


equations forming a discrete-plate model for the steady-state operation. 
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The mathematical developments of this thesis in terms of distilla- 
tion column dynamics or transient models continues in Section M. In 
that section, the discrete-plate steady-state model in Figure 12.3 is 
seen to be a subecase of the discrete~plate dynamic model, and a continu- 
ous-spatial steady~state model en presented) is seen to be a subcase 


of a continuous-spatial dynamic model. 
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CHAPTER 13 


PHILOSOPHY AND DISTILLATION 

This chapter is purely opinion and will be written using the 
first person. lt seems to me that some comments can be made about 
two very general philosopnical desires as applied to the process of 
distillation. The first of these is the seatite to observe and explain 
reality and the second is the desire to achieve profit. Philosophy 
is the study of the principles of reality, and distillation is a 
process for separating components; these two would seem to be almost 
completely unrelated, but, on the contrary, I think a definite and 
crucial relationship exists and propose to present it in this chapter. 
13.1 EXISTENCE EXISTS! 

Tne fundamental axiom of reality is that existence exists, or as 
Aristotle stated, "A is A" (B-28), I exist in this reality as a being 
of voliticnal consciousness with only two choices open to mez: to live 
or to die. I choose to live. I cannot live except by acquiring knowledge 
of one form or another. My only means of acquiring knowledge is through 
my senses. The only tool which I have for the thinking required to 
acquire knowledge from the inputs of my senses is the ability to reason. 
Reason requires the use of logic which is "the art of non-contradictory 
identification" (Atlas Shrugged ~ Ayn Rand), Thus, I wish to acquire 
knowledge about reality by observing it and applying logic to those 
observations. 

One of the many questions I might ask is, “Does there exist a 


defined limit to the amount of knowledge attainable from the observations 
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of my senses and any devices which I may design to aid them?" If 

I wish to observe the operation of a distillation column, can I say 
beforehand that I°1l never be able to understand completely the exact 
mechanisms or motions undergone to achieve the desired separation of 
components? I think that there are several practical resaons why 
the answers to these questions for an individual must be YES, even 
though the philosophical axiom A = A implies that the answers must 
be NO, 

The first practical constraint that I run into is my limited 
physical capacity for knowledge in terms of my limited lifetime for 
acquiring it, But suppose for a moment that I am granted an infinite 
lifespan in which to observe and apply logic to distillation columns. 
I might then attempt to make more specific models of the column based 
upon more accurate observations of its operation. At some point, 
however, I run into another practical constraint: my act of observing 
the operation of the column begins to affect the operation which I'm 
trying to understand. This would be some form of “uncertainty principle" 
applied to distillation column measurements. 

Suppose I don't accept uncertainty; after all, I*°m certain that 
A= A, This supposition means that I can eventually determine a model 
of the operation of a distillation column that tells me to N-— oo decimal 
places exactly what each molecule, atom, or even sub-atomic particle 
is doing within the column at any time, including the effects of my 
observations. Thus, philosophically, the axiom A = A implies that 


unlimited reason is competent to define reality. 
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The certainty of limited life spans of individuals imposes a 
practical constraint on the amount of knowledge about distillation 
attainable by any one man. The assumption that human life will 
propagate forever through time implies that unlimited reason exists. 
Accepting this assumption one can say that man, now and future, has 
the unlimited reason necessary to define reality. What, then, deter- 
mines the portion of the limited facilities of an individual which will 
be applied to the study of distillation? This, then, gets into the 
area of the second philosophical desire, or necessity: to achieve profit 
in order to live, 

13.2 INDIVIDUAL PROFIT MAXIMIZATION 

The amount of effort devoted by any person to the study of distil- 
lation will be in proportion to the benefit (or profit) returned to 
that individual as a result of his efforts, For some this benefit 
may be derived by placing high value upon the self-satisfaction of 
explaining, even partially, the operation of a complex physical system, 
i.e. knowledge for the value of knowledge. Knowledge contains no life 
sustinence, and one who places a very large personal value in knowledge 
must have other means of attaining that profit convertible to life 
sustinence, either by selling that knowledge for food or by devoting 
a portion of his efforts, which could have been spent on studying 
distillation, to the acquisition of food by some other endeavor, 

It is my opinion that each individual attempts to maximize the 
value to him (profit) of the use of his intellect, The individual who 
receives a small (but adequate for life support) income and devotes 


the remainder of his life to watching television is maximizing his 
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value implicitly by placing high value on the endeavor of watching 
television. However, in order for such an individual to have received 
that income necessary for life support, he must have devoted energies 
to some endeavor in which the profits are transformable into food. 
Thus, the conclusion can be made that no matter how an individual 
seeks to maximize the returns from his endeavors some portion of his 
life must be devoted to endeavors which produce something which can be 
transformed into food. 

The question of whether or not there will definitely be people 
throughout time who will devote a portion of their life to studying 
distillation then becomes equivalent to asking whether or not distil- 
lation is an endeavor which produces something which can be transformed 
into food. If distillation is a profitable endeavor, then the relative 
question of how much of an individual's life will be spent studying it 
is answered by defining just how »vrofitable distillation is to the man 
who studies it. 

13.3 PROFIT IN CERTAINTY 

A distillation column has an input, which is a mixture of two con- 
ponents, and two outputs each of which, one is certain, contains a 
higher concentration of a given component of the mixture, It is a 
fact that in the present world the total value of the two outputs is 
greater than the sum of the mixture value and the value given up to 
perform the separation. Thus, distillation IS a profitable endeavor, 
WHY? 

Distillation is profitable because certainty is of very high 


value. The operation of a distillation column produces two outputs 
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of greater certainty than the input. The fact that the composition 
of the outputs is certain is directly transformable into values 
equivalent to food and, as such, is 4 profitable endeavor for an 
individual to use for life support. The existence of profit from 
some endeavor implies that the benefits derived from it are greater 
than the costs incurred. 

What is the cost of certainty? Suppose that I knew a microscopic 
man capable of distinguishing between the two different molecules of 
a binary mixture and capable of deflecting one type of molecule in one 
direction and the other type in another direction. My microscopic 
friend could then be the Maxwell’s Demon of distillation if I let him 
stand in the inlet of the feed pipe to the distillation column and 
bat the lighter molecules upward and let the heavier ones fall, sepa- 
rating the mixture for me. I could then shut down the reboiler and 
condenser and sell the resulting pure outputs with no cost to me in 
terms of energy supplied to the distillation column. In fact, I could 
reason that the oniy effort involved in the entire process is the intel-= 
lectual effort exerted by my smali friend in recognizing which molecules 
are which, Thus, the cost of certainty is the information supplied by 
my friend. 

A supposed fallacy of this argument is often presented as follows. 
The second law of thermodynamics states that entropy is always increas- 
ing in a real, physical process; entropy is then equated to average 
uncertainty, and uncertainty is defined as loss of information. So, 
somewhere down the line, someone loses out; in this case it was, 


according to this argument, my small friend. The Demon's information 
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was being transformed into my profit. The cost of certainty is 
increased uncertainty, if one accepts the concept of average uncer- 
tainty. 

Suppose, once again, that I state that my concept of existence 
denies uncertainty the right to exist for all time. The concept NOT 
is the negation of reality and, as such, cannot exist. No-one will 
ever observe the color Not-Blue because it cannot exist as such. 

The use of the concept of average uncertainty is a denial of the 
statement A = A and can only be rejected as false if A = A is accepted 
as true. The average uncertainty argument implies that uncertainty is 
always increasing; A = A implies that certainty is increasing. The 
two are diametrically opposed. I accept A= A and deny that the 
universe is running down. 

What then is the cost of the separation achieved by my small 
friend if I reject the concept of average uncertainty? The cost is 
the certainty used by my friend in separating the mixture. The profit 
which I derive from the sale of the two separated components results 
from the fact that the two separated components have greater certainty 
and therefore greater value than the value of the certainty used by 
my friend in separating the mixture. Thus, I can pay my friend for 
his services and we'll both be better off as a result of the process, 
Thus, all profit is the result of certainty. 

13.4 CERTAINTY AND KNOWLEDGE 

The only route to certainty is through knowledge. The only route 

to knowledge is through observation and reason utilizing logic. If 


I desire certainty about the operation of a distillation column, then 
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I must turn one on (or read about someone who did), watch it, and apply 
reason to my observations to acquire knowledge about it. My limited 
physical capacity and lifespan dictates that my optimum acquisition of 
knowledge about distillation occurs during the time span up to the 
point where my marginal acquisition of certainty for expenditure of 
concentration begins to become negative. Beyond that point I would be 
better off studying some other profitable venture. Thus, constraints 
must be imposed on the acquisition of knowledge about distillation. 

Tne physical capacity constraints of an individual are realisti- 
cally applied to the study of distillation through judicious approxi- 
mation. I realize that given enough time and intellectual capacity 
I could explain distillation to any desired completeness. However, in 
order to acquire any knowledge about distillation to use in deriving 
profit from my study, without already understanding it completely, I 
must make simplifying approximations and build up a hierarchy of 
knowledge pertaining to the subject. This is exactly what is done, 
and the following sections in this thesis present a theory which over- 
flows with simplifications and approximations. The final test of any 
theory based on approximations must be that the revenue resulting 
from application of the theory is greater than the expenditure of 
effort in developing it. Such a test has yet to be applied to the 
theories presented in this thesis. 
I3.5 WHAT'D HE SAY? 

This chapter has presented my opinions concerning the relation- 


ship between my concept of philosophy and the study of distillation. 
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Tne chapter begins with arguments which position me as a human being 
within the framework of reality as a being who must know to exist, 

The reason I must know is that I must acquire profit (food) to live. 
Distillation is shown to be profitable by arguments which show that 
distillation increases certainty (value) and by arguments which show 
that uncertainty cannot exist if one assumes that existence exists, 
The same arguments are used to invalidate the uncertainty principle 
and the concept of average uncertainty. Finally, the use of simplifi- 
cations and approximations in the study of distillation is justified 


on the basis that a man has limited facilities for use in studying it. 





CHAPTER I4 


A BRIEF REVIEW OF THE LITERATURE OF DISTILLATION 

The volume of literature available pertaining to the study of 
distillation is completely overwhelming! It is easy to see why there 
can be no "Renaissance Man" in modern times, The bibliography pre- 
sented in Chapter Bl of this thesis contains 352 references of which 
202 pertain to distillation column dynamics and control. The 352 
references are minute in number compared to the total literature. 
The 202 references represent a significant percentage of the total 
material available pertaining to column dynamics and control. By 
far the greatest proportion of the total literature of distillation 
deals with column operation, design, and steady-state analysis. 

Distillation is such a general area and has so many inter- 
actions with other areas that it is easy to see why so much can be 
said about it. The areas listed in Chapter B2 represent, in themselves, 
an extensive collection of knowledge, and the list is far from complete, 
In this thesis the term “literature” will refer to this infinitisimal 
subset of 352 references. 

The purpose of this chapter is to present descriptions, comments, 
and notes pertaining to some of the references in Chapter Bl. With 
the exception of some of the theses, each entry in the bibliography 
was very cursorily inspected by the author, and some notes and comments 
were made where deemed worthwhile. This chapter is a connected 
listing of those notes presented in the format of Table B2Z.1. No 


attempt has been made to list complete descriptions, and the author 
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has chosen to use abbreviated language (incomplete sentences, minimum 
articles, etc.) to increase the densite of information with, hopefully, 
minimum loss of comprehension. 

Many of these comments represent this author's opinions and should 
be evaluated as such by the reader. Those references in Chapter B2 
under each area about which no information is given were most likely: 

ae Dynamic response or control theses not read (not available) 
be. Poorly written or invalid (author's opinion) 
c. Not understood by the author after reading, but seemingly 
applicable 
d. Those presenting material covered better in other references 
(author's opinion) 
e. Partially explained by the title 
f. Not really significant to the area but containing pertinent 
material 
g Discussed under another area in which case the area will be 
listed, 
The last name of the first author of each reference is listed for the 
mnemonic convenience of the reader. In the interest of space conserva- 
tion without loss of meaning all area titles and subtitles in Table 
B2.1 have been numbered and brought to the margin. 
T4,.1 GENERAL THHORY OF DISTILLATION 
1. Textbooks 
(B-1) Bennett & Meyers - covers heat, mass, momentum transfer, 
binary and multicomponent distillation, stagewise operations - 


primarily chemical engineering presentation; good presentation 


este 





(G-3) 


(H-10) 


(R=-23) 


(R-21) 


(V-2) 


(A-19) 
(B-15) 


(M-14) 
(0-2) 
(H-2) 


of basic principles of distillation and steady-state analysis. 
Gould = very general process control text; presents detailed 
mathematical methods, models, and analysis for the dynamic 
behavior and control of a large number of chemical processes, 
including packed and plate columns. 

Henley & Staffin - very clear presentation of basic principles 
of distillation. 

Reid & Sherwood = presents in great detail the equilibriun, 
diffusion, thermodynamic, and thermal properties of liquids 

and gases, 

Robinson & Gilliland - detailed presentation of basic principles 
and steady-state analysis. 

Van Winkle - guide to fractional distillation design; considers 
hydraulics, multicomponent aspects, all sites of distillation 
design of plate and packed columns; discusses history of distil- 
lation; see also 3. 

Aris - steady-state distillation analysis. 

Bodman = presents Fortran IV program for optimum design of 
stagewise - ethylbenzene vacuum distillation reactor; economic 
optimization. 

McCabe & Smith - recent text presenting basic principles, 
Oliver - recent text presenting basic eT 

Holland =- uses Thiele-Geddes calculational procedure and 4 
method of convergence; mostly steady-state; "Instead of seeking 
exact analytical solutions for models that roughly approximate 


the actual system, researchers have put a vast amount of effort 
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into the development of iterative procedures in which 
progressively better initial values of the independent variables 


are selected for each successive trial," 


(H-9) ~ 8, 

(H-18) Hengstebeck - Considers many aspects of column design; pre- 
sents basic principles and steady-state analysis; practical 
text, 

(L-26) - 9. 

(C=1) Campbell - steady-state analysis; dynamic analysis using 
signal flow graphs, frequency analysis; process control. 

(P-3) Peters - plant design using detailed plate and bubble cap 
design equations; economics of distillation; optimum reflux 
design. 

(P-10) Pratt - recent text presenting basic principles, steady-state 
analysis, and column design principles. 

(R=22) - 6, 

(R-17) Rosenbrock & Storey - extensive coverage of practical mathemat- 
ical techniques for process dynamics, 

(S-8) Shreve - large number of industrial processes described in 
detail,with diagrams and non-technical descriptions, very little 
mathematics: book represents a lifetime of experience in the 
chemical industry; chemical engineering = unit processes 
(chemical changes) + unit operations (physical changes). 

(S-9)  - 30. 

(S-3) | Sawistowski & Smith - steady-state analysis and calculation. 
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(S-1) - 30. 

(T-1) Treybal - steady-state analysis and column design principles. 

(B-12) Bird - mass transfer in terms of general conservation and 
diffusion equations; some transient analytical analysis. 

(H-5) Henley & Bieber - steady-state example, 3 tray column, q = 2.5, 

benzene - toluene. 

(C-7) Chilton - cost analysis of distillation columns, y = cost 
$ installed per plate or foot, x = diameter squared, then the 
approximate equations for column cost are: 

Plate Columns - stainless - y = 6.47 x0 +63 

Packed Columns - stainless - y = 9,82 ee : 
distillation column maintenance 5-50 $/ft? bubble-plate, 
3-10 $/ft* sieve plate, 5-15 $/ft* for packed columns per 
year in 1960 §$. 

(M-2) = 43, 

(V-3)  VAilbrandt & Dryden - steady-state design of columns; "In 
distillation columns, the pivot point in design is the reflux 
ratio (B in this thesis), which can vary between minimum and 
total reflux. Higher reflux ratios require greater quantities 
of steam and cooling water and a larger column diameter, but 
the column height requirements are lowered. The economic 
reflux ratio is usually 1.1 to 1.2 times the minimum for most 


cases," 


2. kxtensive Bibliographies and Literature Surveys 
Literature Surveys from Industrial and Engineering Chemistry 


Process Control 


(W-16) (W-9) (W-15) (W-8)  (W-19) 
=i 





(W-10) (W-12) (R-25) (W-17) (F-6) 


(W-11) (W-13) (W-18) (G-10) (W-23) 


Distillation 


(B-17) (B-18) (B-19) (F-8) (G-9) (F-7) 


(R-12)-14,; (R-8)-21.; (H-7)-12.; (R-18)-16.; (H-8)-10.; 


(R-17)-1.; (W-26)-41.; (Z-3)-19. 


(G-14) Geddes - gives outline of developments that contribute to 


present knowledge of fractionator design in historical sequence; 
gives comments on present status and future problems; “If 
long-term funds were available for basic research on bubble 
plates, a substantial part of these should be assigned to 


scientific study of the fluid dynamics of plates." 


References of Historical Interest 


3. General Distillation 


(L-19) 


(M-10) 


(L-25) 


(M-9) 
(V-2) 


Lewis - very early (1909) presentation of basic principles; 
uses equilibrium and phase diagrams. 

McCabe & Thiele - the most referenced paper in the literature; 
original (L(25)McCabe-Thiele graphical diagram of steady-state 
column behavior; all previous methods (Sorel's was the first 
in 1899) analytical; algorithm for stepping off tray concen- 
trations on the equilibrium diagram. 

Lewis - early (1922) presentation of stagewise steady-state 
calculations. 

Murphree - explains McCabe-Thiele Diagram in equation forn. 
Van Winkle - brief history of distillation - first recorded 


distillation in Egypt 50 BC, expect earliest unrecorded 


ae 





(U-1) 


(R-29) 


(E-3) 


(C-14) 


(P~12) 


(T-6) 


(G-14) 


distillation 2000 BC, fresh water distilled from sea water 
300 AD, beverage alcohol process first industrial distilla- 
tion process during 11-14th centuries, first books on 
distillation 16th century, stills were differential batch 
type with little reflux up until 19th century, 19th century 
began using steam (1800), bubble caps (1822), continuous 
still (1830), late 19th century first recorded mathematical 
discussions of distillation by Sorel (1899) and Hausbrand 
(1893), then 20th century began mathematics and improved 
distillation (L-19), (M-10), etc. 

Underwood =- the best single presentation of the history of 
distillation in the literature, many pictures and diagrams of 
ancient methods. 

Rodebush - early (1922) plate-by-plate graphical technique; 
preceded McCabe=-Thiele analysis. 

Egloff - review and diagrams of 15th-16th century distillation 
methods and apparatus. 

Cope - early (1932) graphical method using McCabe-Thiele 
stepping type design on the lower part of the equilibrium 
diagram. 

Peters - plate-by-plate calculations on the equilibrium diagram 
(1923) shortly after McCabe-Thiele. 

Thiele & Geddes - referenced many times in the literature, 
graphical and mathematical steady-state analysis and calcula- 
tion method. 

=i 
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(I-1) 


(A-17) 


IBM 705 Program - (1959) steady-state computer program to 
perform Thiele-Geddes calculations; one of the first industrial 
programs, 

Acrivos & Amundson - presents history, basic principles of 
distillation, steady-state analysis using numerical methods 


and eigenvalues. 


4, Dynamic or Transient Analysis 


(M-8) 


(B-29) 


(L-2) 


-7) 


Murphree - (1925) first presentation in the literature of a 
discrete=plate dynamic equation; derives holdup equation, 1 
ordinary differential equation, solves as an exponential. 

Berg & James - (1948) early column transient behavior consid- 
erations; mainly concerned with startup problem, early presenta- 
tion of continuous-spatial model with partial differential 
equations, boundary conditions, and solutions, solves using 
linear equilibrium f(u) = mu + b,.experimentally verified results, 
Lapidus & Amundson - (1950) early transient analysis; general- 
ized the results of (M-1) for countercurrent absorption, showed 
how outlet concentrations could be predicted from the time 
course of the two inlet compositions, assumes f(u) = mu + b 
equilibrium, uses Laplace transforms, poles & zeroes, linearized 
CSE, difference equations, entirely mathematical; calculates 
several transient responses. 

Bartky & Dempster - (1948) early solution to the transient 
system of plate equations, analogous to the classic start-up 
problem except that top reservoir has some holdup as individual 
stages, solution only approximate since compositions derived 


using q # 1.0. 
Se 





(R-2) 


Rose & Johnson - (1953) early development of binary system model; 


relative volatilities, total flow rates, and holdups indepen- 


dent of time, discrete-plate equations solved by Euler predictor 


method, 


(R-12)-14.; (J-1)-21.; (F-5)-19. 


5. Steady State Analysis and McCabe-Thiele Diagrams 


(E-4+) 


(H=-26) 


(C-9) 


(M-18) 


(P-8) 


(S-15) 


(S=20) 


(2-2) 


Eckhart & Rose - steady-state prediction, discrete-plate 
equations, linear equilibrium f(u) = mu + b. 

Hartland - analytical, steady-state comparison, boundary 
conditions. 

Cichelli - steady-state analysis, relates number of plates 

and reflux ratio to the sharpness of separation in binary 
batch distillation, graphs and equations for operation at any 
desired separation are given, many curves, sharpness of separa- 
tion defined in terms of the “pole height," Rayleigh equation 
see 

Mills - List of steady-state computer programs for equilibriun, 
enthalpy, etc., calculations. 

Prausnitz - steady-state computer subroutines for column cal- 
culations for bubble temperature, dew temperature, bubble 
pressure, and dew pressure. 

Surowiec - ideal cascade requires twice the minimum number of 
stages. 

Strand - good development of discrete-plate steady-state 
equations, 

Zuiderweg - general steady-state comparison of device charac- 


teristics and plate efficiencies. 
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(F-3) 


(F-4) 


(H=30) 


(J-7) 


(L-12) 


(A-10) 


(B-26) 


Friday & Smith - discusses mathematically the formulation of 

a solution method for the equilibrium stage steady-state model; 
develops procedure for solving the concentration matrix equa- 
tions which avoids truncation error build up, does not require 
mesh points, works equally well for any number of feeds and 
side streams, and handles nondistributed components. 

Furzer - nonuniform vapor distribution causes reduced efficiency, 
investigates plug flow model and perfectly mixed model, shows 
that maximum reduction in efficiency is halfway between these 
two models. 

Himmelblau - general, recent process mathematical modeling 
text, block diagrams, frequency analysis, matrix analysis. 
Jenson & Jeffreys - steady-state distillation analysis and 
mathematics, numerical methods for solving ordinary and partial 
differential equations, matrix methods, orthogonal function 
theory; primarily a mathematics text. 

Lowenstein - uses normal (Gaussian) probability distribution 
scaled graph paper to greatly increase the accuracy of the 
McCabe-Thiele diagram at the ends, i.e. at very high and very 
low separations. 

Amundson & Pontinen - discrete-plate numerical steady state 
calculations, uses cubic equilibrium relationship, quadratic 
expression for enthalpy, temperature distributions, 4-5 
iterations on 15 plate column; results cited in (M-15). 

Barker - efficiency, design, experimental analysis of bubble- 


cap trays. 


Ito 





(S-12) 


(D-2) 


(R-30) 
(S-24) 


(S-23) 


(S-29) 


(P=9) 


(E-1) 


7) 


(S16) 


sargent - steady-state multicomponent column numerical technique 
and digital computer design; considers individual plate 
efficiencies and equilibrium curve; no dynamics, 

CEP Reprints - collection of steady-state distillation 
papers; articles covering heat and mass transfer, vapor- 
liquid equilibria, packed columns, tray column steady-state 
performance. 

Rose =- experimental evaluation of column steady-state, 
Surowiec - steady-state column design using McCabe-Thiele 
diagrams and discrete-plate equations. 

Stanislas - general steady-state characteristics; plate 
efficiencies, equilibrium, etc. 

Sujata - plate-by-plate, steady-state calculations. 

Floyd = locating feed trays for lowest cost operation, 
steady-state, 

Edmister - true boiling point (TBP) defines distillation 
curve, distillation curve considered continuous and steady- 
state solved by graphical integration technique. 

Treybal = graphical technique for finding plate efficiencies 
based upon McCabe-Thiele diagram, operating lines and equil- 
ibrium curve. 

Srygley - optimum steady-state design in the sense of minimum 
number of plates for desired product purity; uses Thiele- 


Geddes §=-method of convergence, sequential search for optimun. 


(B-1)-1.; (M-10)=-3.; (H-10)-1.; (H-18)-1.; (M-8)-4.3 (M-9)-3.; 


(B-10)-30.; (C-1)-1.; (V-2)-1.3; (G-3)-1.3 (H-5)-1.; (H-2)-1.; 
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(A~17)-3.; (C-14)-3.; (0-2)-1.; (T-6)-3.; (R-5)-16.; (R-21)-1.; 
(T-1)-1.; (S-3)-1.; (H-21)-19.; (I-1)-3.; (L-25)-3.; (A-19)-1.; 
(P-12)-3.; (R-1)-16.; (R-29)-3.; (S-1)-30.; (S-7)-20.; (R-18)-16. 

6. Structural Design 

(L-13) Lowenstein - design of plate sizes using a nomograph and plate 
size "slide rule.” 

(R-22) Rose - minimum cost structural design considerations including 
wind loads, dead weight stresses, longitudinal stresses, 
thickness formulas, numerical examples, drawings of industrial 
columns. 

(J-3) Jones & Van Winkle - experimental analysis of 3 inch perforated 
plate column to determine plate thickness effects on column 
properties, 

(M-11) Manning - structural screens introduced to provide better 
mixing and increase column efficiency. 

(L-12)-5.; (D-11)-16. 

7e conomics and Operations Analysis . 

(M-13) Mitten - economic optimization of distillation by dynamic 
programming, 

(B-15)-1.; (S-8)-1.; (F-9)-5.; (C-7)-1.3 (V=-3)-1. 

8. Thermodynamics 

(H-9) Hougen - Covers,very extensively, thermodynamics applicable 
to distillation. 

(R-23)-1. 

9. Hydrodynamics or Fluid Mechanics 

(P=4)-5,; (V-2)-1.; (B-12)-1.; (R-23)-1. 
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(L-26) Levich - analytical solutions to convection and diffusion 
equations including chemical aspects, 

(H-14) Holm - concludes that vapor flow effects can sometimes cause 
Murphree efficiencies to be greater than unity. 

(F-5)-19.; (T-7)-5. 

(B-4) Bernard - considers sieve tray mixing, foam density, flow 
properties. 

(S-19) Sakata - time for mixing, analytical equipment, mixing pools, 
plug models, and tray efficiencies. 

10. Chemistry 

(B-16) Black - simplified approach to phase equilibria, large number 
of phase equilibria for various compounds presented. 

(H-8) Hala - extensive list of chemical compounds referenced to 
entries in a large bibliography; many phase diagrams. presented, 

feet )-1.; (H-10)-1.; (F-9)-5.; (T-6)-3. 

(H-17) Harper & Moore - experimental paper showing a small still for 
measuring vapor - liquid equilibria lines, several lines given 
as examples. 

(H-20) Howard - concludes that any unsteady-state distillation cal- 
culations should include plate, condenses, and reboiler 
holdups in order to include enough degrees of freedom, 

li. Philosophy 

(A-13) Aris - oftentimes engineer’s rules of thumb are the only tools 
necessary to solve problems adequately, gives examples. 

I4%.2 DISTILLATION COLUMN DYNAMICS 


12. Textbooks 
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(M-1) Marshall & Pigford - beginning of quantitative analysis of 
unsteady-=state operation; many simplifying assumptions used 
to obtain analytical solutions using Laplace transforms, 
assumed equilibrium relationship f(u) = mu + b with m, b 
depending only upon the identity of a component, assumed total 
flow rates and holdups independent of plate number and time; 
considers startup probien. 

(H-7) Holland - uses numerical methods to solve process differen- 
tial equations; uses 6-method for distillation discrete- 
plate equations, example step response of 3-tray column cal- 
culated; good literature survey; book emphasizes numerical 
methods and matrix methods. | 

fea) 1. 

(F-16) Franks - modeling of chemical processes for the purpose of 
control, plate equations, partial differential equations. 

13. Theses 

(M-15) Mohr - considers each of two sections of a column independently, 
determines step response for each input stream using IBM 705 
and discrete-plate equation, each of the response curves is 
then fitted with a two-time constant exponential expression, 
these expressions manipulated into general transfer function 
using Laplace transform variable p, time constants and gain 
factors of the column are then correlated with steady-state 
parameters from McCabe-Thiele diagram; analysis based on linear 


equilibrium f(u) = mu + b; assumes binary mixtures, constant gq, 
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(R-6) 


(B-32) 


constant fluid rates, total condensation, negligible holdup, 
constant liquid holdup each plate, reboiler holdup equivalent 
to that on each plate, q =1, perfect mixing, unity plate 
efficiencies; found that major time constant was independent 
of the type of disturbance, gains and secondary time constants 
vary with disturbance type and with initial steady-state, 
predicted responses of large columns more sensitive to changes 
in values of selection parameters than small columns; column 
time constants increase with H/L, degree of separation, degree 
of nonlinearity of equilibrium curve, number of plates. 
Romagnoli - hyorid simulation of discrete-plate equations; 
butadiene distillation plant, 2 towers of 49 trays each of 
bubble type, constant pressure, constant holdups, water 
ignored, vapor holdup ignored, constant Q, En = 0.7, Francis- 
weir formula, micro assembly language and PDP=-1 used; con- 
clusion was that hybrid setup worked faster than purely digital 
computation. 

Brosilow & Tanner - develops ard analyzes methods for computing 
and optimizing countercurrent staged processes using distil- 
lation as an example; formulates models in terms of both 
discrete and continuous spatial equations; solves continuous 
models using economic cost criterion using gradient and 
Lagrangian non-linear programming methods; found equations 
much easier to solve when models were not restricted to 
integral values of stages, ailowing more general optimization 


methods to be used such as Lees’ method; two-point boundary 
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(G-5) 


value problem solved by imbedding it into an initial value 
problem in time; found that when Lees’ algorithm was modified 
to generate tridiagonal matrices the rate of convergence 
increased substantially. 

Gaydos ~- develops generalized digital computer program and 
model for a single bubble-cap plate; model includes multi- 
components, hydraulics, non-ideal vapor-liquid equilibriun, 
heat transfer to and from each stage, provisions for feed and 
take off streams; computation time limitations were encountered 


in the simulation. 


14. Reviews, Bibliographies, and Literature Surveys 


(A-3) 


(R-12) 


Archer & Rothfus - presents survey of the 1955-1960 dynamic 
behavior literature; discrete-plate equations developed, 
startup and transition between steady states, batch operation 
of plate and packed columns, and process control are all 
discussed in terms of their respective literature; frequency 
analysis. 

Rosenbrock - surveys the history and present developments of 
discrete and continuous distillation and heat exchanges 


models; best presentation of this material in the literature. 


(H-7)-12.; (R-8)-21.; (Z-3)-19.; (L-3)-15.; (T-2)-16. 


(S-33) 


Lehigh Symposium - discrete-plate models, transient response, 
and control of distillation columns; extensive bibliography of 
149 references on column steady-state, dynamic, and control 


aspects; frequency response analysis. 


Additional surveys, see 2. 
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Dynamic Models or Solutions 


Discrete Plate Equations 


15. 


Frequency Analysis or Laplace Transform Solution 


(A-3)-14,; (R-12)-14.; (H-7)-12.; (R-8)-21. 
(Z-3)-19.; (S-33)-14.; (T-2)-16. 


(L-3) 


Lamb, Pigford & Rippin - oscillations in tray compositions 
resulting from input oscillations in either feed composition 

or reflux flow are calculated using analog computer for 16- 
tray column; found frequency response at low frequencies like 
simple first-order process and at high frequencies having 
interference patterns and large phase lags; equations linearized 
about steady-state, E=1, frequencies from 0.001 to 1 radian/ 


tray holdup time, 5-tray column frequency response also obtained. 


Transient or Time Analysis 


Numerical Solution 


16. Digital Computer 


(D-H) 


Distefano, May, & Huckaba - (1967), discrete-plate dynamic 
model solved for a sequence of upsets in which the next step 
occurs before the transients of the previous one have died out; 
solves large system of equations by Adams-Moulton-Shell (AMOS) 
finite difference technique on IBM 709, computation time was 
12 mine, expect nybrid steup would eed this up, 12-plate 
column, spacing 1 ft, 10" diameter, samples of every 3rd 

plate, 7 runs for different types of steps and pulses; 

oriented toward prediction for feedforward control; transient 


times 40-80 min., experimental data agreed with computer 
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calculations to within 5%, error blamed on truncation in 


predictor-corrector methods. 


(B-32)-13.; (R-4)-4.; (R-20)-31. 


(D-8) 


i-3) 


(L-5) 


(D-15) 


Distefano - numerically solves discrete=-plate equations by a 
large number of different methods and then compares then, 
(D-10) presents stability aspects. 

Huckaba - numerical solution by IBM 650 of transient response 
to binary 1Z-plate column, experimentally confirmed; uses 
nonlinear discrete-plate equations, inputs are step changes in 
heat input to the reboiler, feed composition, and simultaneous 
changes in feed composition and reflux ratio; computation time 
5 min. using fourth-order Runge-Kutta for starting and fifth- 
order modified Adams for continuing. 

Luyben = uses set of linear perturbation type differential 
equations for characterizing dynamic behavior; experimental 
results verified for acetone-benzene system, 1.8 ¢q < 2.2; types 
of disturbances were changes in feed composition, feed rate, 
top tray reflux rate, bottom tray vapor rate, equations 

solved by analog computer, 

Duffin & Gamer =- defines model for multicomponent distillation 
which includes secondary effects due to column hydraulics, 
holdup, and delay effects in boundary conditions; errors in 
generated system response at low values of elapsed time 
usually due to inadequate models; used general discrete- 

Plate model and numerical integration scheme, conclude that 


model is valid. 


-5/- 





(R=16) 


(R-9) 


(R-10) 


Rosenbrock = surveys available solution methods for the discrete- 
plate equations, decides to use digital computer for speed and 
generality; describes numerical method and program, computa- 
tion time 5 min. for 5 plates up to 100 min. for 300 plates. 
Rosenbrock - develops discrete plate equations, discusses 
possible methods of solution, almost identical to (R-16). 
Rosenbrock = discusses relative advantages of two computer 
programs for solving the equations in (R-9); first program 
has f(u) fed in as table of 101 values, linear interpolation, 
forward integration used, first program finally rejected due 
to limitations; first program required solution of large 
system of simultaneous equations, second program developed to 
eliminate this by solving step-by-step by equating slopes; 


routine included to evaluate df(u)/du at discrete-points. 


(R-11) Rosenbrock - discusses the accuracy of computer programs in 


(R-10), extends application to multicomponent systems; used 
equilibrium curve as f(u) = u + 0.02, E = 1.0; concludes that 
the most promising method for calculating transient response 


is digital computer. 


(A-4)=40.; (R-17)-1.; (S-33)-14.; (S-16)-5.; (S-20)-5. 


(M-3) 


Mah, Michaelson, & Sargent - dynamic behaviour of multistage 
systems described by large sets of non-linear first-order 
differential equations; discrete-plate equations then linearized 
but shown to be inconsistent due to linearization, step-by- 

step procedure using exponential function is proposed, 


numerical integration technique, reviews all standard numerical 
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(T-2) 


(W-1) 


(P-2) 


procedures, uses tridiagonal matrix, finds eigenvalues, shows 
that use of linear equilibrium f(u) = Ku also leads to physical 
inconsistency, 

Tetlow, Groves, & Holland - develops a generalized model 

which accounts for the effects of channeling, transfer lag, 
mixing, and mass transfer in unsteady-state, multicomponent, 
discrete-plate distillation; model tested on a large number 

of numerical examples; considers plug flow and perfect mixer 
holdups, large tridiagonal matrix results, @-method of con- 
vergence used; use of the generalized model in the analysis 
of control systems is discussed. 

Wood & Armstrong - derives a linearized model of the discrete- 
Plate equations using f(u) = mu + b equilibrium; solves for 
step response of feed composition; comparison with experimental 
results shows that the model is only valid for moderate values 
of time after the step and cannot be used as the column 
approaches final steady state. 

Peiser & Grover - presents a model similar to (H-3) including 
the effects of heat and mass transfer and tray hydraulics; 
predicted that unsteady-state prediction can be used to solve 
several significant problems in multicomponent distillation 
which are not evident from steady-state analysis; numerical 
computations carried out on a digital computer simulating 


column open and closed loop control. 
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(D-9) 


(L-18) 


(D-6) 


(S-6) 


iv-9) 


7-1) 


(G-4) 


(R-5) 


(D-11) 


Davison ~- solves large systems of x = Ax + Bu equations; 
calculates the poles and zeroes of the system, then solve 

for a few of the more significant variables in terms of poles 
and zeroes. 

Luyben - discrete-plate equations solved for transient 
response for use in feedforward control. 

Davison ~- discrete~plate equations solved using matrix methods 
for the transient response of a column due to pressure varia- 
trens for use in control. 

Sargent ~- discrete~plate equations solved numerically using 
matrix methods, 

Thorogood ~- discrete-plate equations solved using Runge-Kutta 
methods. 

Yesberg & Johnson - demonstrates use of a resistance network 
analog to solve the absorber problem of (A~1) and (L-2), 
discrete-plate equations linearized, time derivative represented 
by a backward finite difference to produce a set of simul- 
taneous algebraic equations which are solved on an IBM 650 by 
matrix inversion, 

Greenstadt ~ discrete-plate equations solved by Newton's 
method, 

Rose, Sweeney & Schrodt - solves startup problem for ternary 
mixture using discrete-plate equations, Lewis-Matheson method 
used. 

DiLiddo & Walsh - pulse column considered as a series of 
stages, 3 plate ideal model, 9 plate real model, numerical 


solution on IBM 605. ee 











(R-1) 


(R-19) 


(L-27) 


(w-4) 


(R-18) 


Rose, Johnson & Williams - (1950) discrete-plate equations 
solved by Euler predictor method; model assumes relative 
volatilities, flow rates, and holdups independent of time. 
Rose, Johnson & Williams - (1951) similar to (R-1), early 
papers showing pictures of IBM cards, plate equations solved 
by numerical methods. 

Lowe ~ discrete stage equations solved on digital computer, no 
distillation. 

Waggoner & Holland - discrete-plate equations solved using 
Simpson's rule, 3 point corrector to approximate the integrals 
in the component material balances, results presented in table 
form. 

Rose & Williams - (1950) early plate-by-plate solution of 
discrete model, similar to (R-19) (R-1), shows wiring of 


computer control panels. 


17. Hybrid Computer 


(F-1) Franks - solves several countercurrent problems, several of 
these were for distillation. 

(F-13) Frank & Lapidus - discrete-plate equations solved by hybrid 
computer using one integrator for each plate. 

(R-6) 13.3 

(F-12) Frank & Lapidus - hybrid simulation particularly useful for 


nonlinear partial differential equations, uses repetitive 


analog operation, memory, and first-order lags. 
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18, Analog Computer 


(B-23) 


e-1) 


(R-3) 


(G-12) 


(L~5)-16. 


Bowman & Clark - uses linear equilibrium f(u) = ku which was 
found to be valid near the top of the column where volatile 
component concentration is small, discrete equations for 20 
and 30 plate columns wired up on analog computer; when column 
was Switched from total reflux to finite reflux ratio, an 
almost instantaneous drop in overhead composition occurred, 
from this point on the overhead was independent of the time 

on total reflux and dependent only on the stillpot composi-~ 
tion at the end of the total reflux period; this suggests that 
column can be divided into two independent sections, first 
would be period on total reflux, second would be behavior at 
finite reflux; linear equilibrium does not provide exact rep- 
resentation of real column, provides guide to the degree and 
direction of change within the column. 

Pigford, Tepe & Garrahan - solves unsteady-state equations for 


batch distillation of a binary mixture using a mechanical 


analog computer called the “differential analyzer"; assumes 


constant relative volatilities, vapor & liquid flow rates, 
Rose & Williams - demonstrated use of an analog computer to 
obtain transient response and deSign a controller for a 5=-plate 
column, uses analog computer with Pade delay circuits; large 
number of problems solved to determine the best controller for 
maintaining the composition of the distillate constant under 
composition and thermal variations of the feed. 

Grover & Peiser ~ analog computer solves plate equations for 


control, stability aspects considered. 
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19. 


(D-3) 


(G-2) 


Analytical Analysis 


Davidson - uses mechanical analog and Rayleigh's method 

to determine eigenvalues for finite series representation of 
plate column transient behavior; plate-type stripping column 
fed at top, no bottoms takeoff, assumes linear equilibriun; 
equilibrium concentration on any plate proportional to 

exp (-8T), where g depends upon q and number of plates N; 
Rayleigh method used to give approximate 8; solves example 
model of Taylor Diffusion type, found that first term in 
series was most important. 

Gilliland & Mohr - analytical analysis of discrete-plate 
equations using two-time constant exponential model of (M-15); 
digital computer solves for transient response, then the two 
time constants are determined from the results and used to 
develop transfer functions; by use of transfer functions, 
the responses of several columns to step changes in feed 
composition were predicted and compared with the responses 


calculated by numerically solving the discrete-plate equations, 


(A-3)-14.; (M-8)-4.; (M-1)-12.; (R-8)-21.; (S-31)-43. 


(R-7) 


Rosenbrock - funcamental “disturbance trapping" paper; char- 
acterizes the departures of a distillation column from its 
steady-state by a quantity D which measures the rates of change 
of composition on all the plates and increases whenever the 


column is disturbed from its steady-state; result is generalized 
D=25 | a(H,u_)/at| 
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(R-32) 


(W-2) 


to multicomponent systems; eee. become “trapped: in 

the column at points where equilibrium curve slope equals 
operating line slope, rates of change of composition do not 
decrease as would be expected, numerical computation of behavior 
in these cases is difficult; presents an electrical analog 

to the discrete-plate equations. 

Rosenbrock < discrete-plate equations, matrix methods, energy 
considerations. 

Wilkinson & Armstrong - considered response of a column at 
total reflux to a change in bottom vapor composition, used 
binary mixture, linear equilibrium curve; response predicted 
by analytical analysis was in good agreement with that observed 


experimentally. 


(F-16)-12.; (G-3)-1.; (H-20)-10.; (B-7)-4.; (B-12)-1. 


(B-27) Balasubramanian - analytically solves one differential equation 


(C-2) 


(2-3) 


5) 


for a one-plate still. 

Cullinan - considers transient start-up problem using matrix 
methods for analytical solution, uses f(u) = mu + b. 

Zykov - analytical solution of discrete-plate equations for 
transient analysis of multicomponent column; good bibliography 
of Russian literature. 

Foss, Gerster & Pigford - assumptions of complete mixing or 

no mixing lead to inaccuracies in distillation, paper attempts 
to establish the nature and extent of mixing and to develop 
calculational methods to account for its effect on plate 
efficiency; mixing experimentally determined by use of 


tracers and measurements of residence times; simplified calicu- 
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(H-21) 


lation procedure presented which affords a rapid means of 
computing plate efficiency under all mixing conditions. 
Hassett - solves for transient behavior using equilibriun, 


steady-state, and McCabe-Thiele diagran,. 


Continuous Spatial Equations 
20. Frequency Analysis or Laplace Transform Solution 


(J-6) 


(M=27) 


(H-4) 


7) 


(W-5) 


(D-1) 


Jafri, Glinski & Wood - continuous system transient response 
with control using time constants and transfer function analysis. 
Majumdar - solves CSE using Laplace transforms and lineariza- 
tion, applies to Clusius column. 

Hoerner & Shiesser - frequency and time responses using Laplace 
transforms, linearized CSE, linear equilibrium; gives model of 
Taylor Diffusion type, 

Sellers & Augood - transients in a liquid hydrogen separator, 
uses Laplace transforms, exponential characterizations, rate 

of approach method. 

Ward = uses Laplace transforms and frequency analysis to cal- 
culate the time behavior of dynamical systems. 

Douglas & Rippin - uses linearized equations and sinusoidal 


inputs, considers system in terms of chemical oscillators, 


21. Transient or Time Analysis 


(J-1) 


Jackson & Pigford - (1956) digital computer solution of linear- 
ized CSH model for startup problem, plots of composition 
throughout column as a function of reduced time are presented; 


linearized CSE and linear equilibrium f(u) = Au gives equation 
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of Taylor Diffusion model type; took 3 hours of computation 


time on 1BM 701 to solve for transient curves on several trays. 


(H-4)-20.; (D-3)-19.; (M-1)-12.; (C-2)-19.; (L-26)-9, 


(K-2) 


Koppel - solves heat exchanger/chemical reactor equation of a 
form possibly similar to CSE; 


du=[1+r(t)]du+ Pfr t+ or(t)j us, 
gt Ox 


(O-1) Osborne = linear equilibrium f(u) = mu + b used, equations 


(R-8) 


(S-13) 


considered continuous in time and differenced in theoretical 
stage; concludes that the real cause of numerical instability 
problems is in theoretical stage direction differencing; 
numerical solution obtained with very large theoretical stage 
step, one theoretical stage, and very small time steps, this 
took care of instability problems; Fortran IV computer program 
for max. 6 components, max, 60 trays described. 

Rosenbrock = calculates transient responses; describes labor 
needed to solve CS# models; describes control aspects and 
theory; presents good bibliography. 

Stone & Brian - detailed numerical methods described for 
solving CSE type equations; solves Taylor Diffusion model 
equation as an example; CSE is a subcase of equation (2) 


which is a general form of convective transport equation; 


O Cee) du] ~o_ [V(x,t,u)f(u) | =9u 
ax gx gx at 


several different types of numerical methods discussed and 


compared; analysis derived applies only to linear equations; 


Bre 





(M-28 ) 


(K-3) 


(2-7) 


(R-14) 


recent results indicate that the desirable features of 
solutions obtained by the new equations for linear problems 
are to a large degree found in solutions of nonlinear problems. 
Montroll & Newell =- exact solution of nonlinear differential 
equations which describe the time dependent behavior of multi- 
stage cascade separating processes of two very similar non- 
linear species, Rayleigh separation law postulated for each 
stage; linearization of nonlinear second-order portion differ- 
ential equations is discussed; analytical expressions in terms 
of exponentials and eigenvalues are developed; equilibrium 
curve is approximated by f(u) = u + cu(l - u). 

Kermode & Stevens - several nonlinear continuous models 

solved on an analog computer. 

Powers - numerical solution to Taylor Diffusion type partial 
differential equation with two-point boundary conditions and 
one initial condition. 

Ruckenstein - analytical analysis of convective diffusion 
(Taylor Diffusion) type equations; extensive analysis of 


applicable transformations and linearization methods. 


(B-29)-4.; (G-3)-1. 


-5) 


(H-25) 


Tsang - solves heat equation using eigenfunction expansion 
(orthogonal), evaluates lst 3 values numerically; other 
analytical solutions presented also, Bessel functions, 

Herron & Van Rosenberg - uses a mesh and centered difference 
method to numerically solve convective transport equations with 


two-point boundary conditions. 
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(J-2) 


(J-5) 


(L-20) 


(B-20) 


(W-14) 


(B-25) 


(C-8) 


Jury - solves heat equation using analog computer, uses 
memory, solves repetitively. 

Jackson - numerical solution and optimization of partial 
differential equation models. 

Lapidus - transient response and control of chemical reactors 
using continuous models. 

Bedingfield & Drew - heat and mass transfer expressed by the 
Same equations. 

Woodle - analogy between distillation and heat transfer, some 
equations. 

Brian - transient response using a continuous, Taylor 
Diffusion type model. 

Crank - diffusion in different geometrics with different 
boundary conditions, finite difference methods, diffusion 
and chemical reaction; use transformations for equations witn 


variable diffusion coefficients. 


Experimental Transient Behavior 
22. Frequency Response 


(H-15) 


(H-23) 


(A-6) 


Henley - frequency response techniques for experimental 
analysis of transient behavior. 

Hutchinson & Shelton - frequency response techniques using 
correlation functions. 

Armstrong & Wilkinson - carried out experimental work to verify 
theoretical computations and methods of (R-16) and (R-9); 


studied behavior of 21 plate CH), ~ CCl), separator subject to 


B63= 





(H-24) 


(W-28) 


step changes in feed and reflux on column composition; long 
time agreement was better than short time after disturbance; 
results summarized by transfer functions having the form of 
pure time delay followed by a linear log, both time constants 
are functions of plate number. 

Haagensen - experimental results derived from frequency 
response using matrix techniques, 

Woods = experimentally determined controller settings for a 


continuous system operating dynamically. 


23. Time Response 
(H-3)-16.; (L-5)-16.; (R-30)-5. 


=>) 


(B-3) 


ei-15) 


Baber & Gerster - determines experimental transient response 
of column to changes in liquid and vapor rate, demonstrate 

the applicability of discrete-plate equations for predicting 
measured response; responses to step inputs presented in 
tables and graphs; equations solved by analog computer; 
results confirmed validity of model; linear perturbation types 
of equations predict satisfactorily the transient behavior. 
Baber - similar to (B-5) but earlier; experimental response 

of 5-tray, 2 ft. bubble-cap column; tests made over a range of 
gas and liquid rates nearly up to the flooding point and for 
tower pressures up to 5 atmospheres; average difference between 
predicted and experimental results was 13%, indicating that 
Simple perturbation equations are valid models, 

Rademaker = tested an ethylene-ethane splitting column to 


provide data for checking a general theory of column dynamics; 
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data summarized in graphs and experimental accuracy discussed. 

(D-12) Davies - "hidden transients" can have significant effects on 
tray efficiencies, 

24h, Cyclic Distillation 

(A-14) Atkeson - one of the first papers (1957) to indicate that 
cycling or nonequilibrium operation can increase mass transfer 
rate greatly; attempts to explain physically. 

B2.3 DISTILLATION COLUMN CONTROL 

25. Textbooks 

(A-12)-40.; (G-3)-1.;(C-1)-1. 

(B-13) Buckley - very practical, industrially oriented process control 
text. 

(K-5) Koppel - matrix theory, optimal control, sampled-data control. 

(A-11) Athans & Falb - mathematical theory of optimal control; many 
examples and problems; very clear presentation; extensive 
bibliography. 

26. Theses 

(B-2) Beecher - presents method of dynamic control system synthesis 
using calculus of variations; provides extensive Fortran program; 
example is 100 plate butadiene/butiene-2 column. 

(B-32)-13. 

(G-1) Gordon-Clark - uses matrix theory to adjust dynamic response 
of process to some required response before applying conven-~ 
tional control; uses linearized model of a 5=plate binary 
column, control and measure composition on each plate; main 


computation difficulty is finding eigenvalues of matrix; 
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Fortran program, form matrix, form into tridiagonal matrix, 
obtain complex roots of real polynomial by Bairstowe's method, 
iteration using quadratic factor, obtain roots of quadratic 
factors, return; major drawback is large number of variables 
required. 

27. Extensive Bibliographies and Literature Surveys 

(A-3)-14.; (R-8)-21.; (R-12)-14.; (S-33)-14. 

Conventional Control Systems 


28. Digital Control 
(H-29) Hanson, Duffin & Somerville - provides extensive Fortran 


computer programs for control. 

(B-6) Buster - closed loop control applied to oil fractionating 
system, required large memory and clever programming to 
complete calculations in real time. 


29. Hybrid Control 


see 17. 

30. Analog Control and Instrumentation 

(H-16) Haines - large number of diagrams of possible column control 
configurations with explanations. 

(G-3)-1.; (B-13)-25.; (C-1)-1.; (G-12)-18. 

(B-22) Buckley - basic column control strategy, diagrams of control 
loops, linearized models. 

(B-10) Bauer & Orr - uses McCabe-Thiele diagram to derive best 
operating iines for Bontrons 

(B-31) Boyd - mostly nonmathematical and non-diagramatical discussion. 

(P-6) Pink - describes control using analog computer, no mathematics. 
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(S-1) 


(L-4) 


(C-12) 


(T=-11) 


(N-29) 


(P=14) 
(R-13) 


(S-9) 


(H-1) 


Shinsky - process controls based upon time response analysis; 
distillation columns hard to control because (1) many trayed 
towers slow to respond to control action, (2) separation 
requires many variables, (3) on-line analysis not aiwaye 
available, (4) distillation units are usually last in the chain 
of processing operations, and (5) factors affecting separation 
not readily interpreted in terms of control system requirements. 
Lupfer & Parsons - describes a control system designed to reduce 
the effects of changes in flow rate and feed composition on 
column operation, uses predictive or feedforward control, 
dependent outputs of the process are controlled by measuring 
one or more inputs (which generally cannot be controlled), 

and then the controlled parameters are changed as required 

to achieve the desired output. 

Ceagiske - comparison of transient and frequency response 
methods for control of linear chemical process systems. 

Tivy - control discussion, no mathematics. 

Moczek - discusses effect of transient behavior and dead time 
on control, very little mathematics, 

Phillips - control descriptions, no mathematics. 

Rijnsdorp - discusses feedforward & feedback control with 
distillation column as an example. 

Strobel - describes in detail theory of optional and electrical 
measuring devices; great deal of theory presented, photometers, 
spectral analysis, wave theory, etc. 

Harriott - transient and frequency analysis used for column 


control. 
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Control Systems Using Dynamic Models 

31. Digital Computer Control 

(A-12)-40.; (D-4)-16.; (L-18)-16.; (D-6)-16. 

(R=20) Rosenbrock - application of automatic control theory to 
chemical processes has not led to the same improvement of 
performance as it has in the control of mechanical and elec- 
trical systems, paper attempts to explain why; complexity of 
chemical processes, lack of suitable measuring equipments, | 
remote objectives i.e. out of 1000 variables only 10 are of 
interest, modes not always easily separated; proposes matrix 
method in which important modes are picked out for control. 

(S-33)-14. 

(C-11) Cadman, Rothfus & Kermode - uses matrix methods and frequency 
response techniques for design of multicomponent feedforward 
control system. 

32. Hybrid Computer Control 

(D-13) Dahlin & Nelson - uses hybrid computer and matrix maximum 
principle for optimal control. 

33. Analog Computer Control 

(J-6)-20.; (K-2)-21.; (L-3)-15.; (R-3)-18.; (H-15)-22. 

(L-14) Lupfer & Oglesby - elaborate analog controller described, 
instrumentation and feedback schemes. 

(L-6)  Luyben & Gerster - studies effectiveness of feedforward 
control for 10 and 40 tray columns, performance of overhead 
and bottoms controller determined by analog simulation and by 


experimental tests on a 10 tray - 2 ft. column; concludes that 
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relatively simple feedforward controllers appear adequate for 
distillation; for small input disturbance, a linear model can 
be used to determine controller transfer functions. 

(W-27) Williams & Harnett - uses frequency response analysis, plate 
equations, first order lags; describes various control schemes. 

34. Optimal Control 

(L-9) Lapidus - nonlinear optimal control, quadratic performance 
criteria, Riccati equation and solution, etc. 

(B-11) Brosilow & Handley - optimal control of the overhead composi- 
tion of a distillation column, integral squared error criterion 
on disturbances; experimental analysis on 5 inch column with 
15 trays and 3 bubble caps per tray} control system behaved 
well in spite of model inaccuracies, 

(B-32)-13.; (A-11)-25.; (D-13)-32.; (K-5)-25.3 (S-33)-14.; 

(J-5)-21. 

35. Distributed or Modal Control 

(D-9)-16.; (J-6)-33.3 (S-31)-43.; (F-16)-12. 

(G-13) Gavalas - eigenvalue solutions for distributed parameter 
steady state. 

B2.4 MATHEMATICS AND COMPUTATION 

36. Ordinary Differential Equation Theory 

(H-11) Hartman - $20.00 text, ordinary differential equations in 
all aspects from a pure mathematics standpoint, theorem-proof 
presentation. 

(B-14) Birkhoff & Rota - very good presentation of transformations and 
eigenfunction expansions for two-point boundary-value problems 


and Sturm-Liouville problems. 
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(1-2) 


(H-12) 


Ince - if one were ever restricted to only one look on ordinary 
differential equations, this would have to be the one. 
Hildebrand ~- very practical applied mathematics text, useful 


theory of infinite series expansions to solve DE's, 


37. Partial Differential Equation Theory 


(G-7) 


(F-11) 


(B-21) 


Garabedian - theory of lst and 2nd order PDE's, basics of 
integral equation theory, very few examples or problems. 
Forsythe & Wasow ~- practical numerical methods text, useful 
for solving PD's. 

Berg & McGregor ~- very clear presentation of basic principles 


and solution techniques for PDE's, 


(O-1)-21.; (F-12)-17.; (F~-16)-12. 


(G-6) 


(W-7) 
(H~28) 


Gurel & Lapidus - extensive discussion of stability of ODE's 
and PDE‘s. 

Webster ~- large number of practical example-problems solved. 

Hildebrand ~- newest Hildebrand text, seems to be as practical 
and understandable as the previous texts; numerical methods 


for solving PDE's. 


38. Integral Equation Theory 


3) 


(P-4) 


(L-8) 


Tricomi - best text on integral equation theory in literature, 
complete, readable, very few examples. 

Petrovskii - good presentation of separable kernel theory and 
use of algebraic equations to approximate integral; several 
classifying examples and problems. 

Lovitt - very clear presentation of basic theory with large 


number of easily worked and instructive examples and problems. 


~7 5- 





(H-13) 


(D-5) 


(M-6) 


(S-11) 


(v-H) 


Hildebrand - very explicit presentation of methods to convert 
from IE's to DE"s with examples. 

Davis - large number of methods and equations solved which 
are not found anywhere else in the literature. 

Mikhlin - very readable text, terminology sometimes different 
from standard U.S., coveres about same material as (L-8). 
Smithies - rigorous text, presents theory of IE's in terms 


of Lebesgue integration and L, spaces; Riemann integration 


fi 


and R,, Spaces are subcases of Lospaces. 


Volterra - original classic in theory of integral equations. 


(G-7)-37.; (W-7)-37. 


(W-6) 


Whittaker & Watson - very clear, but abbreviated, presentation 
of IE theory; extensive presentation of infinite series 


expansions. 


39. Mathematical Transformations 


(T-4) 


(2-1) 


Tranter - basic tneory and application of different types of 
integral transforms. 
zemanian = recent text, extensive theory of transformations 


and integral transforms. 


(S-31)-43.; (M-16)-44.; (R-14)-21.; (C-8)-21. 


4O. Matrix Mathematics 


(A-4) 


Acrivos & Amundson - transient solution to column discrete 


AL solution; this paper refer- 


equations by matrix metnods, e 
enced very often in the literature; paper originally brought 
the subject of matrices to the attention of chemical engineers, 


wide variety of chemical engineering problems solved in this 
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(A-12) 


paper by matrix methods including the absorber problems of 
(A-1) and (L-2) a 
Amundson - good presentation on eigenfunction expansions; 


matrix theory. 


(H-13)-38.; (J-7)-5. 


41. Numerical Solution Techniques 


(B~30) 


(R=244) 


Berry & DePrima - develops iterative procedure for the 
determination of the eigenvalues and eigenfunctions associated 
with the solution of Sturm-Liouville problems in a finite 
interval; presents and discusses convergency of an iterative 
scheme different from “sweeping” or Rayleigh-Ritz; presents 
numerical example. 

Ralston = very practical and widely referenced text, many of 
the programs in the Scientific Subroutine Package (I-3) use 


methods of this text. 


(D-8)-16.; (H-28)-37.; (F-11)-27.; (B-32)-13. 


(F-2) 


(K-1) 


(L-10) 


Fox - somewhat dated but clear presentation, verified with 
many numerical examples, of methods for finding boundary- 
value solutions and eigenvalues. 

Kenneth & McGill - gradient methods, general numerical methods, 
some theorems on existence and uniqueness of boundary-value 
solutions. 

Liu =- numerical solution by finite-difference method; shown 

to be explicit, stable, more accurate than Crank-Nicholson; 


solves heat equation and several nonlinear examples. 





(M-22) McGinnis - solves Taylor Diffusion type equations using 
Runge-Kutta techniques; describes numerical methods for BV 
problems. 

(S-13)-21.; (C-8)-21.; (A-10)-5.; (P-8)-5. 

(H-19) Hamming - very practical and extensive numerical methods text. 

(M-4) McCracken & Dorn - application of numerical methods, flow 
charting, basic Fortran, programming principles, ODE's and 
PDE’s numerical methods, etc. 

(R-4) Rose, Johnson & Williams - used both analog and digital 
computer to solve plate equations for a 7 plate binary column; 
found that the time required for the column to change from one 
steady-state to another after an abrupt change in feed composi- 
tion is a strong function of the magnitude of change, q, reflux 
ratio, N, feed tray location; numerical methods used described 
briefly. 

(R-29)-3.; (P-12)-3.3 (H-25)-21.; (J-7)-5.3 (A-17)-3. 

(L-11) Lee - invariant imbedding approach; classical methods use 
quasi-linearization; invariant imbedding considers a family 
of problems from zero to the duration of the original problen, 
by imoedding, solve for the missing conditions for 2-point 
BV problem; solves Taylor Diffusion type equation as an example. 

(D=-9)-16.; (M-2)-43.; (R-17)-1.; (J-5)-21. 

(N-1) Naylor - numerical solution techniques. 

(W-26) White - method for numerically calculating eigenvalues and 
eigenvectors of large dimension matrices; good bibliography. 


(D-15)-16.; (S-12)-5. 
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42. Boundary Value Problems 
(B-8) Boyce & DiPrima - basic text on practical solution techniques 


and theory of boundary and eigenvalue problems. 

(B-9) Beltrami & Wohlers - very theoretical presentation; existence 
and uniqueness theorems for boundary value problems. 

(V-1) Villadsen & Stewart - new collocation methods given for 
solving symmetrical boundary-value problems using orthogonality 
conditions to select collocation points; accuracy is shown to 
be comparable to least squares or variational methods, 
calculations are much simpler; applications given to one- 
dimensional eigenvalue problems and to parabolic and elliptic 
PD's; collocation methods are special techniques for solving 
integral equations numerically. 

(L-11)-41.; (H-11)-36.3; (K-1)-41.; (B-14)-36.; (M-22)-41.; 

(B+32)-13.; (F-2)-41.; (H-12)-36.; (H-25)-21. 

43. Higen-Values, Vectors, and Functions 

(S-31) Singer - state variable transformations and matrix methods 
to select significant modes and eigenvalues of multivariable 
systems. 

(M-2)  Mickley, Sherwood & Reed - numerical solutions of PDE's 
using finite differences applied to stagewise processes; 
orthogonal functions and infinite series solutions of PDE's, 

(B-30)-41.; (A-12)-40.; (G-13)-35.3 (H+11)-36.; 

(H-13)-38.; (I-2)-36.; (B-14)-36.; (D-3)-19.; 

(S-13)-21.; (V-1)-42.; (W-26)-41.; (B-5)-22.; 


(H-12)-36.; (J-7)-5. 
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Wi, Special Functions 


(A-18) 


(L-23) 


(M-16) 


(R=28 ) 


Abramowitz & Stegun - complete and extensive analytical, 
numerical, and graphical presentations of special functions. 
Lebedev - practical user's text showing mathematical properties 
of special functions and examples of usage. 

Magnus, Oberhettinger & Soni - states and proves many of 

the mathematical properties of special functions, 


Rainville - similar to (L-23). 


45, Computation and Computer Programming 


(0-H) 


(I-3) 


Organick - basic principles of Fortran IV in textbook form 
with large number of examples. 

Scientific Subroutine Package - set of over 250 subroutines 
for performing standard numerical manipulations such as 
matrix inversion, integration, differentiation, expansion in 
functions, least squares curve fitting, roots of polynomials, 
etc.; (I-3) is in Fortran IV but recently (1969) a reduced 


package in PL/I has been xeleased. 


(M-4)-41.; (P-8)-5.3; (H-29)-28.; (M-18)-5.; (B-15)-1. 
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SECTION 2 


MODELS OF BINARY DISTILLATION COLUMNS (M) 


Ml THE DISCRETE-PLATE EQUATIONS (DPE) 
M2 THE CONTINUOUS-SPATIAL EQUATION (CSE) 
M3 SOLUTION TECHNIQUES FOR THE CSE 


M4 LINEAR APPROXIMATIONS TO THE CSE 


“I HAVE HARDLY EVER KNOWN A MATHEMATICIAN WHO WAS CAPABLE 
OF REASONING." — PLATO 


“PLATO WAS A FOOL!" — JACKIz GLEASON 


THIS SECTION PRESENTS DISCRETE AND CONTINUOUS MODELS FOR THE 
COMPOSITION BEHAVIOR OF A BINARY PLATE DISTILLATION COLUMN. THE 
CONTINUOUS-SPATIAL EQUATION (CSE) IS DEVELOPED FROM THE DISCRETE- 
PLATE EQUATIONS (DPE) BY CONSIDERING THE PLATE NUMBER TO BE A 
CONTINUOUS VARIABLE. SOLUTION METHODS FOR THE CSE ARE PRESENTED, 
INCLUDING TRANSFORMATIONS, GENERAL SOLUTION TECHNIQUES, NONLINEAR 
APPROXIMATIONS, AND LINEAR APPROXIMATIONS. AS A FINAL RESULT THE 
LINEAR POLYNOMIAL-COEFFICIENT MODEL (LPCM), WHICH IS THE SUBJECT OF 
SECTION 3, IS DEVELOPED AS A LINEAR APPROXIMATION TO THE CSE, 
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CHAPTER M1 


THE DISCRETE-PLATE EQUATIONS (DPS) 


The development of the discrete dynamic model of a binary plate 
distillation column is the subject of this chapter. The development 
begins by defining the symbols to be used in the derivation. Then 
the individual plate equations will be developed by considering com- 
ponent mass balances on a typical plate, Finally several comments are 
made referencing the literature pertaining to the discrete-plate 
equations. The developments of this chapter are a simplified version 
of those presented in reference (G-3). 

The symbols to be used in the derivation of the individual plate 

equations are presented in Table Mi.l. These symbols are to be 
applied to the characteristics of the idealized, typical plate of the 
lower or stripping section of the column as shown in Figure M.l. 
The derivation of the equation for the rectification or upper section 
are identical except that the upper flow constants must be used. The 
author has attempted to keep these symbols consistent with those used 
in Chapters Il and I2 and especially in Table 12.1. 

The dynamic model is a system of equations developed for each 
plate by utilizing the basic principle of conservation of mass stated 


in Equation Ml.1. The four possible mass balances which can be applied 





to the plate of Figure Ml.l1 are listed below. Three of these mass 
balances, numbers 2, 3, and 4 are satisfied by the assumptions stated 
in Table Ml.1 for the individual variables, These assumptions restrict 


~82~ 





=u) 


Oo DB RR 


= 


Concentration of the lighter component in n-th plate liquid 
(lbm mole lighter component/lbm liquid) 


Concentration of the lighter component in the n-th plate vapor 
(1bm mole lighter component/lbm vapor) 


Relative volatility (dimensionless); assumed constant 

Feed rate (lbm liquid/hour); assumed constant 

Distillate rate (lbm liquid/hour); assumed constant 

Bottoms withdrawal rate (lbm liquid/hour); assumed constant 
Liquid holdup on the plate (lbm liquid); assumed constant 
Vapor holdup above the plate (lbm vapor); assumed constant 


Liquid rate in the rectification section (upper section) 
(lbm liquid/hour); assumed constant 


Liquid rate in the stripping section (lower section) 
(lbm liquid/hour); assumed constant 


Vapor rate in the column (lbm vapor/hour); assumed constant 
Upper reflux ratio, L,/V, (lbm liquid/lbm vapor) 

Lower reflux ratio, L,/V, (1bm liquid/lbm vapor) 

Feed plate index (integer) 

Internal plate index (integer), len<N 


Total number of plates in the column (integer) 


Murphree plate 


fficiency (dimensionless), defined by 
Equation I2.1, 1 


assumed for all plates 
Portion of the feed which adds to the lower liquid rate; 


L = Ly + qgF ; further developments use q = 1, saturated 
itquid“feed 


Table Mi.1 - LIST OF SYMBOLS USED IN DISTILLATION 
COLUMN DYNAMIC MODELS 
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1. Liquid phase mass balance of the lighter component, 

2. Liquid phase overall mass balance. 

3. Vapor phase mass balance of the lighter component. 

4, Vapor phase overall mass balance. 
the overall validity of the model greatly but are utilized to simplify 
the developments which follow in later chapters, It is eventually 
expected that the generality of the Linear Polynomial-Coefficient Model 
(LPCM) presented in Section 3 can be utilized to generate solutions 
which are nearly as accurate in representing the essential nature of 
the column transient behavior as would be a completely general model 


using four equations per plate, 
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Migure Ml,1 ~ TYPICAL RECTIFICATION PLATE WITH SYMBOLS 











The mass balance of the lighter component in the liquid phase 
gives equation Ml.2 for the upper section and Ml,3 for the lower 
section. The steady-state portions of these equations are seen to 
be equivalent to the discrete-plate steady-state equations in 
Figure 12,3. 


H du, 
at 


V (f(y,.,) - f(u,)] + el RR M1.2 


H du, . 
av 7 Y (fCu_y) - f(u,)) + fy -) 2.3 

The overall binary distillation column dynamic model utilizing 
the discrete-plate equations can now be presented as Figure Ml.2, 
where the effects of the condenser and reboiler have not been included, 
It can be seen that this model consists of N-nonlinear ordinary differ- 
ential equations with boundary conditions. This particular model is a 
highly simplified version of the more general models (using all 4 
equations with varying holdups and flow rates) usually considered in 
the literature (See B2,2 ~ Discrete Plate Equations). The solution 
techniques to be applied to the continuous-spatial equation developed 
in later chapters can also be applied to more general models; but for 
the purposes of presenting the LPCM and its soliution techniques, the 
model of Figure M1.2 is sufficient. Wnen accurate quantitative behavior 
is desired, the discrete-plate model is almost always the model 
employed in the control system applications described in the literature 
(B2.3 = Digital Computer Control). 

The reader is referred to the very extensive literature in 
Chapters Bl and B2 for further information pertaining to discrete- 


Plate models and solutions. This chapter has presented a very brief 
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development of the discrete-plate model (DPM) for the purpose of 
utilizing it in the development of the continuous-spatial model (CSE) 
which is subsequently to be linearized to the linear polynomial- 
coefficient model (LPCM) for which an analytical solution procedure 
is developed, The next chapter presents the development of a contin- 


uous-spatial model (CSE) from this discrete-plate model, 


Model Equations 


Rectification Section |e a A 0 = 9) 
H du, 
ae Y PQ) - PQ) + B, ay] 
Feed Tray n=k 
sl Gol 


“k - 
Fal [f(u,3) - flu) ] + byway Lju, + Fug 


Stripping Section ge Mar 1 
H du, . 
a7 Y CQ) - PQQ] + Day] 


Equilibrium 


f(u,) = an 


ee (a-i)u, 
End Conditions 


F(u_) i Top n=N 


I 
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Feed and Flow Conditions 
L, =L,+aF or B, = B, + aF/V 
V=L +0D 
u 
W= Ly - V 


Figure M1.2 = THE DISCRET@-PLAT™ ZQUATIONS FOR_A BINARY 
PLATE DISTILLATION COLUMN 
Se 








CHAPTER M2 


THE CONTINUOUS-SPATIAL EQUATION 


This chapter presents the steps in the procedure used to transform 
the discrete=plate equations of Figure Ml.2 to the continuous-spatial 
equation (CSE) which is a continuous dynamic model of the concentration 
behavior of a binary plate distillation column. In general this con- 
version represents a significant step away from the quantitative accuracy 
of the discrete models because of the series expansions and approxima= 
tions involved, Usually, however, the solutions of the discrete models 
require comparatively large amounts of time to calculate, often on the 
order of half (about 15 minutes) of the major concentration time con- 
‘stant (B2.2 = Digital Computer). This makes such models of limited 
usefulness in control systems which utilize the model to predict 
transient response in order to correct for transients ahead of time. 

For these reasons a model is sought which requires less computation 
time to predict the dynamic response of the distillation column, 

The continuous-spatial equation (CSE) may be such a model. The 
remainder of this thesis will be devoted to developing, investigating, 
approximating, and solving models of the continuous~spatial type. These 
models are not expected to behave quantitatively as near to the actual 
column response as do the discrete models, but it is anticipated that 
some sacrifice in numerical accuracy will result in a great savings in 
computation time and therefore better prediction for use in control. 

The literature pertaining to continuous models is somewhat limited. 
The best single presentation on this subject in the literature of 
Chapter Bl is by Rosenbrock (R-12), It would seem that the most general 
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form of continuous model solved analytically in the literature is the 
Taylor Diffusion Model (Chapter L4) which results from approximations 
which, for the case of distillation, have been shown not to represent 
accurately the behavior of the column (G-3) (R-7). The most general 
form of continuous model solved numerically in the literature is the 
polynomial-class nonlinear model (Chapter M3). This model is solved 
by numerical methods which treat the continuous variable as a set of 
discrete points, and therefore the solution takes nearly as much time 
as solving the discrete equations (B2.2 ~ Transient or Time Analysis). 
It can be expected that most of the time savings resulting from solving 
continuous models will be a direct result of extensive analytical 
analysis and judicious approximation. This will be the objective of 
the following chapters in this thesis. 

Tne development of the continuous-spatial equation begins with 
the representation of the plate number n by a continuous-spatial 


variable x given in equation M2.1. The variable x is continuous from 


x = en M2.1 


Where: n = plate number 
g = plate spacing (ft.) 


ttom to top of the column. The general discrete-plate equation is 


represented here from Chapter Ml as equation M2.2. The transformation 
H du, = vf(u,_3) - Vf(u,) + Lui, 7 Lu, M2.2 
dt 
Where: 2 ee u(t) 


from plate number to spatial variable means that Uy, is redefined as 


BB 





in equation M2.3 and M2.4, Using these functions, equation M2.2 now 


becomes the partial differential equation of M2.5. 
u(t) =x, U) M2.3 
U41 (+) = u(xtg,t) M2.4 


H @u(x,t) = vr [u(x-g,t)] - vf [u(x,t)] M2.5 


i 
+ Lu(xtg,t) - Lu(x,t) 


The continuous-spatial equation M2.5 is as exact as the discrete- 
plate equation M2.2, only it is in a different form. Equation M2.5 is 
too general for any analytical solution, thus, the terms in x-g and 
xte are each expanded in a Taylor Series (G-3), and up to the second 
order terms are retained as in equations M2.6 and M2.7. If these are 


now substituted into equation M2.5 the result is given by equation M2.8. 
vf [u(x-g,t)] = ve [u(x,t)] - ev@_f [ulx,t)] 
ox 


2 O2 
, a oe 1 Rea ee M2.6 


Lu(xtg,t) = Lu(x,t) + eb Of rc, C) | 
ax 


+ g°L 9? u(x,t) M2.7 
2 9x° 


Hdu = g2 9? [ilut vf(u)] + so [Lu - Vf(u) ] 
gt 2 Ox? 9x OO 


Equation M2.8 can be written in a somewhat neater form with the 


use of several variable changes. Dividing through M2.8 by V 
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and considering t in all of the above equations to be t = ty» where 
Ee is the column time variable, then the transformation to a new time 


variable t is given by M2.9. Similarly the continuous spatial variable 


of the column is considered to be X,, and the change of variable on 


x is given by equation M2.10. If the condenser and the reboiler are 
x= x /6 M2,10 


now considered as the O'th and N + 1*st plates respectively, then the 
ranges of the variables can be scaled as in M2.11 and M2.12, with the 


column height represented by C, (ft.) in M2.13 for equally spaced plates. 


Oe x. =< CU Meee 

aac = Th 
trax =< 1.0 M2.12 
= (N+l)g Mees 


If all of these operations are applied to equation M2.8, the result 
is equation M2.14, 


dusl 92 (Bu+ f(u)] +92 _[Bu-f(u)] 2.14 
a 2 Ox? 4x 


The complete model of thé distillation column composition dynamic 
behavior using continuous spatial equations (CSE) can now be presented 
as Figure M2.1. This chapter has presented the steps in the procedure 
used to convert the DPE model of Chapter Ml to the CSE model of Figure 
M2.1. Tne CSE model is a highly nonlinear partial differential equa- 


tion of second order with nonlinear boundary conditions and as such 
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requires that approximations be made before any analytical results 
can be expected, Linear approximations to this model are the main 
subjects of the remainder of this thesis. The next chapter describes 
possible solution techniques applicable to the CSE and develops the 


Linear Polynomial-Coefficient Model (LPCM). 


Model Equation 
du=107 [But f(u)] +9_ [Bu - f(u)] 
Btu 2 ax” ox 


Where: B= L/V Upper Xp < x < 1.0 
Eee &/V Lower 0 < xX < Xp 


Xp feed tray location 
Equilibrium 


f(u) = son 


Boundary Conditions 


Qu = f(u) -u 


B Qu -Of(u) =F (u-u,) 
ax ox V 


For: B= L/V and B, = L,/V 


Ba 


ox 


=f(u) -u 


Heewcemic. | — THh CONTINUOUS-SPATIAL EQUATION (CSE) 


DYNAMIC MODEL 
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CHAPTER M3 


SOLUTION TECHNIQUES FOR THE CONTINUOUS-SPATIAL EQUATION (CSE) 


This chapter presents a discussion of possible solution techniques 
applicable to the CSE, with particular emphasis on a nonlinear ee 
imation to the LSE. The basic equation of the CSE is here repeated 
for convenience as equation M3.1. The CSE is a second-order, nonlinear, 
partial-differential boundary-value problem with nonlinear boundary 
conditions. An outline of possibie techniques for solving the CSE 
is presented in Figure M3.1. In this chapter Parts A, B, C, and E 
of Figure M3.1 are discussed very briefly, and Part D is examined with 
the aim of developing the Linear Polynomial-Coefficient Model (LPCM), 


which is examined in detail in Section 3(1). 














M31 





Pi 


at 


5 [Bu + f(u)] + s [Bu = f(u)] 


As stated previously, any distillation column model would consist 
of two separate CSE's of the form M3.1: one for the rectification or 
upper section of the column and one for the stripping or lower section. 
The two separate CSE's are connected by the boundary condition at the 
feed tray, which represents the dividing point between the lower and 
upper sections of the column. Any complete solution for the dynamic 
composition behavior of a distillation column would therefore involve 
combining the solutions from two CSE's and their appropriate boundary 
conditions. The central problem to be considered now is that of 


solving the CSH, 


Sone 





In order for the reader to get some idea of how complicated the 
CSE really is, the equation is presented in expanded form as equation 
M3.2. Considering Part E of Figure M3.1, the general theory of 
partial differential equations (B2.4 - Partial Differential Equation 
Theory) provides little help in finding the solution to the CSE. 
Tne general theory presents extensive and valuable results pertaining 
to existence of solutions, characteristic analysis, and methods of 
solution for linear, semi-linear, and quasi-linear equations, none of 


which can be applied to the completely nonlinear CSK,. 


du = ee Of(u)] d2u + (1 2?#(u) , au 
Z gue 9x C2 Qu Ox 


at M3.2 
- Of(u) + B] Qu 
gu 9x 
Where: u = u(x,t) , B= constant = L/V 


io 
‘h u) = j- 


Concerning Part A of Figure M3.1, it would seem to be equally 
marlrcuLlt to Simulate the entire set of discrete=-plate equations as to 
simulate the CSE by any of the three methods listed. In fact, since 
the CSE is a partial differential equation, use of digital, hybrid, 
or analog equipment to simulate it requires that the CSE be broken up 
into discrete-spatial parts. Some time savings might result in appli- 
cation of an analog computer for the nonlinear spatial portion evaluated 


at discrete times, 
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A. Computer Simulation of the CSE 
1. Digital Simulation 
2 Hybrid Simulation 
3. Analog Simulation 
B, Application of Transformations to the CSE 
1. General Transform Theory 
ae Frequency Transforms; Laplace, Fourier, etc. 
be Integral Transforms 
2. Change of Variables in the CSE to produce: 
ae More Easily Solved Nonlinear Equation 
b. Linear Equation 
c. More kasily Approximated Equation 
C. Application of Some Results in Pure Mathematics to the CSE 
1. Contraction Mappings and The Fixed Point Theorem 
2. Techniques for Finding the Fixed Point 
D. Approximation Techniques for the CSE 
1. Noniinear Polynomial Approximations 


2e Linearization and Linear Polynomial-Coefficient 
Models (LPCM) of the CSE 


BE. Application of the Theory of Partial Differential 
Equations to the CSE 


Figure M3,.1 - POSSIBLE SOLUTION TECHNIQUES FOR THE CSE MODEL 


Solution of the CSE by application of some form of transforma- 
tion to the variables, Part B in Figure M3.1, is one area which might 


offer rewards with further investigation. Of course, frequency 


Sone 





transforms such as Laplace or Fourier can be eliminated from consider- 
ation immediately in the case of the general CSE because they depend 
upon linearity properties. However, these techniques might be effec- 
tively applied to the linearized CSE, and in fact, they are used for 
analytical and numerical solutions in the literature (B2.2 - Continuous 
Spatial Equations = Frequency Analysis.) 

In the application of a general transformation, an equation of the 
form M3.3 is sought which will make the CSE a more easily solved non- 
linear equation, a linear equation, or a more easily approximated 
nonlinear equation. In addition, a variable transformation of the form 


F(s) = r f(u)K(s,u)du M3.3 


O 


of equation M3.4 could be applied to the CSE, The resulting complete 
set of possible transformations on dependent variables is then given 


by M3.5. Transformations of this type are often used successfully for 
u = h(v) M3.4 


F(s) = r f[h(v) | K(s,v)dv M3. 5 


solving problems in fluid mechanics and heat transfer. A search of 
these areas for transformations applicable to the CSE or some nonlinear 
approximation to the CSE might prove rewarding. 

part C of Figure M3.1 presents an area of analysis which may 
possibly apply to a theoretical examination of the CSE, The major 
consideration here is the use of the Fixed Point Theorem of Brower 


(S-21) (H-11) to prove such properties as existence, uniqueness, and 
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continuity of solutions to the CSE, A description of the Fixed Point 
Theorem depends upon the concept of a contraction mapping. The appli- 
cation of a contraction mapping to the CSE depends upon considering 
the CSE as a nonlinear differential equation of the form of equation 
M3.6 at each instant of time. 

du = f(x,u) M3.6 

dx 

The next step in applying the Fixed Point Theorem is the conver- 

sion of M3.6 to an integral equation (B2.4 - Integral Equation Theory) 


presented in equation M3.7. Then the right hand side of M3.7 can be 


toc) = u(x.) “ in f[z,u(z) jdz MB era 


considered as a general functional transformation (mapping of 
functions) A(v) given by equation M3.8. Equation M3.7 can next be 
written as a functional mapping of the function v onto the function 


u as in M3.9. 
A(v) = va ie f(z,v)dz M3.8 
O 
A(v) = u M329 


Now, any functional transformation A(u) which satisfies equation 
M3.10 for any two functions u, (x) and uy (x) is defined to be a contrac- 


tion mapping. The Fixed Point Theorem of Brower (1912) then guarantees 


| A(u, ) ~ A(uy) jem U, ~ Uy M3.10 
Nuebos 0 auMac 1 


for any contraction mapping that there exists a function u, the fixed 
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point in the mapping, such that M3.11 is satisfied. In other words 


the function transforms to itself. Equation M3.11 is equivalent to 
A(u) =u M3.11 


equation M3.7, and thus the Fixed Point Theorem proves existence of 
a solution to this special case of the CSE if the CSE when written as 
a transformation can be shown to be a contraction mapping. 

The answer to the question of whether the CSE represents a con=- 
traction mapping or not could only be determined by further detailed 
investigations of the CSE. However, since in the original column the 
variable u represented the concentration of the lighter component in 
the liquid, it could certainly be expected that u(x,t) would always 
remain in the range of M3.12. Any distillation column model which 
has negative or greater-than-unity concentration solutions cannot be 
valid. In addition, all of the analytical solutions to approximate 


versions of the CSE in Chapter L4 satisfy equation M3.12. 
Ofeu(x,t) < 1,0 M3.12 


The Fixed Point Theorem offers valuable information concerning 
the existence of a solution to a general nonlinear differential equa- 
tion, but it says nothing about the techniques for finding that solu- 
tion. It is expected that general solution techniques for an equation 
involving transformations and contraction mappings would have to be 
found in the mathematics literature (general topology) and applied 
to the CSE as a special case. It is highly unlikely that any such 


techniques could be applied to the CSE without some simplifying 
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assumptions or approximations being made. 

Thus, the analysis of the CSE now proceeds to Part D of Figure 
M3.1, beginning with nonlinear polynomial approximations to the CSE, 
The essence of nonlinear approximations lies in the assumption that 
the equations generated by approximating the coefficients of equation 
M3.2 using the initial steady-state distribution will be valid for the 
transient behavior of the column. The validity of this assumption 
can only be determined by solving specific examples by numerical 
methods and comparing them to the CSE numerical solutions, 

One possible nonlinear approximation technique would be the use of 
n-th degree polynomials to represent the coefficients of equation M3.2, 
The steady-state equation would then be a polynomial class (D-5, Ch.8) 
nonlinear ordinary differential equation of the form of equation M3,13. 
The resulting nonlinear approximate model of the CSE would be given 
by equation M3.14, where the Ps and q. are determined by approximating 
the coefficient functions of u and the equilibrium relationship by 


polynomials of sufficient degree for desired accuracy. 


A(u) d?u + B(u) du + C(u) duy? + Diu) = 0 M3.13 
dx 


adx® dx 
Mm a 
Wheres A(u) = A,u 
1=0 
n 
B(u) = x Bu 
1=0 
i 
C(u) = C.u 
i=0 
q st 
Bir e= >. Du 
(u) ee) 
Ass By» Cys D, are functions of x 


298% 





¢] 


u= 9? u u : 
gu oe [Py (u)] + e [P,(u) ] M3614 
Where: P, (u) = a pw 
Po(u) = = q,u 


The nonlinear approximated version of the CSE presented in M3.14 
is an analytically unsolvable equation in general. That M3.14 can 
be solved using less computation time than that necessary for the CSE 
is very doubtful. It may be that M3.14 for specific cases might be 
more easily programmed on an analog computer, but further investiga- 
tion for specific cases would be required to determine so. 

This chapter has presented a brief discussion of a variety of 
different methods which could be applied to the CSE, The results of 
the ee scicet ton are as expected; in order to proceed with any 
analytical analysis of the CSE, the equation must be linearized. 

The development of linearization techniques is presented in the next 


chapter (M4). 
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CHAPTER M4 


LINEAR APPROXIMATIONS TO THE CONTINUOUS-SPATIAL EQUATION (CSE) 


This chapter presents the steps in the linearization procedure 
to convert the CSE to a linear partial differential equation with 
spatially varying coefficients. The spatially varying coefficients are 
then approximated by n-th order polynomials in x, forming the Linear 
Polynomial-Coefficient Model (LPCM). The most important characteristic 
of linearized models is that the variable u(x,t) in the linearized 
model represents the incremental distribution resulting from the 
linearization about steady-state operation. 

Tne CSE is written in equation M4.1 in terms of Une) where 
u(x,t) $s the column concentration of the lighter component in the 
liquid. If M4.1 is linearized about the initial steady-state operation 
using equation M4.2, then the resulting linearized CSE model for u (x,t) 
representing the column deviation from the initial steady-state distribu- 
tion resulting from transient boundary condition variations is given by 
M4.3. The details of the linearization procedure are presented in 
Figure M4.1. Equation M4.3 can then be written as a linear partial 

Qu, 197 (Bu, + f(u,)] + [Bug-f(ug)] M41 

ar ax ox 

u(x,t) = u, (x) + u(x,t) M4, 2 
differential equation of second-order with spatially varying coeffi- 


cients by expanding the terms as in equation M4.4, 


% 


u = 1 33 [Bu + m(x)u | a O- [ Bu - m(x)u | M4. 3 
Z ax 


> 
| 
> 
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Development of the Linearized CSE 


u) =190? u uJl+a u- £(u)]) - du 
A(u) +e + f(u) | oe f(u) | - 








A(u,) = A(u,) +@a | (u) 
ou fu, 
| 
ese = Us 
A{u,) = A(ug) = 0 
Ox | (u) = 0 
ou 
OA =10? [(B + m(x))u} +2_[(B - m(x))u] - du = 0 
gu fu, 29x? 9x ae 
Development of the Boundary Conditions 
A. (u) = (a ) f(u) + u = 0 ne ede (0) Sifoyo 
ox 
Expanding similarly to above: 
(ae ae) (m(x)u) + u = 0 
ox 
[am(x) - m(x) + 1] u + m(x) Ju = 0 x = 1,0 Top 
dx ox 
A(u) = BO_[u - f(u)] - F (u - up) = 0 X= Xp, Feed 
ox V 
Expanding as above: 
-[B dm( x) ct FI u + [1 = n(x) ] du = ur(t) 
ax V ax V 
x = xX Feed 
Where: up(t) = up (t) - up; Upper’ and Lower 
Ap(u) = (B, 2 +1) u - f(u) x = 0,0 Bottom 
Xx 
Expanding as above: 
[1 - n(x)]J u+ B, Qu = 0 x = 0,0 Bottom 
Xx 


Figure M4.1 - DETAILS OF THE CSE LINEARIZATION 
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ar 2 re] dx dx 


+ [1 d?m(x) - dm(x)] u MY, dy 
2 dx? ax 


du = 1 | m(x) ] d2u - m(x m( x u 
1 LB + m(x)]e%u + [B (x) +4 du 


It is very important to realize the nature of the function m(x). 


The function m(x) is defined in equation M4.5 and represents the slope 


gu 
u = u, (x) 

of the equilibrium curve as a function of x. In equation M4.3 the 
constant B represents the slope of the operating lines on the McCabe-~- 
Thiele diagram of Figure 12,2 in the upper and lower regions. The 
function m(x) represents the slope of the equilibrium curve in those 
same regions. Almost all approximation techniques applicable to 
equation M4.3 use these characteristics and assume that m(x) is given 
by an equation of the form of M4.6 or M4.8 (R-12) (G-3) (B2.2 - 
Analytical Analysis). With the assumption of equation M4.6 equation 


M4.3 becomes equation M4.7 which is the Heat Equation. 
m(x) = constant = B M4, 6 


du 
gt x” 

It is rather obvious that assumption M4.6 reduces the linearized 
equation M4,.3 to the simplest diffusion equation, thereby neglecting 
most of the behavior which is unique tc distillation. Distillation 


is a diffusion process, and it is somewhat reassuring to find that 


the model M4,3 reduces to a diffusion model. However, it is certain 
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that a model of the form of M4.7 could sot be used to accurately 
predict column behavior for the purpose of controlling the column, 
Equation M4.7 is solved analytically and numerically for B = 1.25 in 
Chapter L4 and represents the "bottom of the ladder" in any hierarchy 
of approximate models of the CSE, 

The next step upward in complexity results from assumption M4.8 
resulting in equation M4.3 becoming M4.9. Equation M4.9 is the Taylor 
Diffusion model equation and is discussed in detail in Chapter L4. 
This equation is still too approximate to use for any control applica- 
tions and so, a more sophisticated representation of equation M4.3 is 
sought. 


m(x) = constant # B Mm4.8 


du = Pig?u + Pau 
at Ox? ax M4.9 





Where s Ps (B + m)/2 


Po = (B - m) 


Suppose the spatially varying coefficient functions in equation 
M4.3 are approximated by n-th degree polynomials in x. ‘Then the 
resulting complete model would be the Linear Polynomial-Coefficient 
Model (LPCM) presented in this chapter as Figure M4.2 and in 
Chapter Ll as Figure Ll.1. The boundary function constants As; in 
Figure Ll.1 are equivalent to the coefficient functions, evaluated at 
the boundaries, in Figure M4.2. 

The central purpose of this thesis is to suggest that the LPCM of 
Figure M4.2 can be used as a distillation column dynamic model of 
sufficient accuracy and rapid solution as to be useful in controlling 


NOS 


a 





the column and to present an analytical solution technique for it. 


This is the subject of Section 3(L). 


Model Equations 
Qu =9?_[P)(x)u] +9_ [P,(x)u] 
gt gx gx 


Where: P, (x) - _P4X” = [B + m(x) |/2 


n al 
P, (x) = y 44% = [B - m(x) ] 
Upper and Lower Equations, One for Each Section 
of the Column 


Boundary Conditions 
[an (x) - m(x) +1] u+ m(x) du =0 x = 1,08 lop 
x 


x 


-[B an(x) + F] u+ (2 ~ m(x)] Ou = =F u(t) 
V Ox V 


ax 


x= xX Feed 
(Upper and Lower) 


[2 - m(x)]u+ B,du = 0 x = 0,0 Bottom 


ox 


Lgure M4.2 - THE COMPLETE LINKAR POLYNOMIAL- 
COEFFICIENT M 
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SECTION 3 


THE LINEAR POLYNOMIAL - COEFFICIENT MODEL (L) 


Ll ANALYTICAL SOLUTION OF THE LPCM BY INTEGRAL EQUATION TECHNIQUES 
L2 THE PARTIAL LINEAR POLYNOMIAL BOUNDARY VALUE PROGRAM (PLPBV) 

L3 DETERMINATION OF THE LPCM FOR EXAMPLE DISTILLATION COLUMNS 

L4 ANALYTICAL SOLUTIONS TO APPROXIMATED COLUMN EQUATIONS 


L5 OTHER SOLUTION TECHNIQUES APPLICABLE TO THE LPCM 


"DAZZLED BY THEIR ABILITY TO DO ELEMENTARY THINGS AT TREMENDOUS 
SPEEDS AND TO PUT THESE TOGETHER IN STRUCTURES OF DAUNTING COM- 
PLEXITY, SOME HAVE ALLOWED THE TERM ‘GIANT BRAINS® TO GAIN 
CURRENCY AND, SEDUCED BY THE SIREN SONG OF SO SENSELESS A 
SOBRIQUET, HAVE SURRENDERED THEIR BIRTHRIGHT OF RATIONAL THOUGHT 


FOR A POTTAGE OF PUNCHED CARDS." - R, ARIS (A-13) 


THE CENTRAL PURPOSE OF THIS THESIS IS THE DEVELOPMENT, 
PRESENTATION, SUGGESTED SOLUTION TECHNIQUE, AND EVALUATION OF 
THE LINEAR POLYNOMIAL - COEFFICIENT MODEL (LPCM) OF THE DYNAMIC 
BEHAVIOR OF A BINARY PLATE DISTILLATION COLUMN, THE FUNCTION 
OF THIS SECTION IS THE PRESENTATION AND UTILIZATION OF THIS MODEL 
AND THE DEVELOPMENT OF AN INTEGRAL EQUATION SOLUTION TECHNIQUE 


FOR IT. 
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CHAPTER L1 


ANALYTICAL SOLUTION OF THE LPCM BY INTEGRAL EQUATION TECHNIQUES 
The purpose of this chapter is to present a somewhat involved 


analytical solution of the Linear Polynomial-Coefficient Model (LPCM) 
defined by equation Ll.1 with boundary conditions L1.2, L1.3, and 
L1.4. The LPCM is a linear, second-order, parabolic, partial- 
differential, two-point, boundary~value problem and, as such, cannot 
be solved analytically or numerically with complete generality (B-8). 
The best that can be done analytically is to develop a solution 
technique which allows reasonable assumptions and approximations to 
greatly simplify the analytical manipulations leading to a simplified 
solution. The best that can be done numerically is to develop and 
employ programs specificaily related to specialized cases of the LPCM 
and then to utilize those programs on digital, hybrid, or analog 
equipment to compute a simplified solution. In either case, the 
solution may or may not truely represent the behavior of the distilla- 
tion column from which the model was developed. 

Equations of the form of the LPCM are usually solved numerically 
by a variety of different techniques. Numerical solutions to the 
LPCM are described in Chapter L5 and a survey of the literature in 
the Bibliography, Chapter Bl, pertinent to solving the LPCM is 
presented there. No attempts have been made to solve the LPCM in a 
purely numerical fashion in this thesis, although the suggestion can 
be made that such attempts might prove to be exceedingly profitable 
in terms of greatly reduced solution computation time when using the 


LPCM as part of a column control system. 
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The solution technique to be described in this chapter is a 
sequence of analytical manipulations which lead to an expression 
which is then to be evaluated numerically. The technique begins with 
expressing the response in terms of steady state and transient 
portions. The technique of separation of variables is then employed 
for the transient partial differential equation. The spatial 
differential equation portion resulting from the separation of 
variables is then solved by converting it to a Liouville Normal- 

Form equation and then to a homogeneous Fredholm II-integral equation. 
The total solution is then expressed in terms of the eigenvalues and 
eigenfunctions resulting from the solution of the integral equation. 
The anticipated advantage occuring from these analytical manipulations 
is that the final numerical evaluation in this technique may take 

much less computation time than solving the LPCM purely numerically. 
Llel TRANSIENT AND STEADY-STATE EQUATIONS FROM THE LPCM 

Tne Linear Polynomial-Coefficient Model is defined by equation Ll.1l 
with boundary conditions L1.2, L1.3, and L1.4. This model represents 
the deviation from steady state of the concentration of the lighter 
component in the distillation column subject to step changes in the 
end-point compositions. The complete presentation of the development 
of this model has been given in Section 2(M). 

The first step to be taken in the analytical solution of Ll.1 
is the realization that, at any given instant of time, there are 
three separate spatial concentration distributions implied by the 
model, First of all, there is the initial steady state distribution 


w, (x) in the column prior to the step inputs. Since the LPCM 
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£ LP, (x)u] + 0_[P,(x)u] oa du 
x ox ot 





Wheres u = u(x,t) 





i 
Py(x) 2 px 





n 4 
P(x) #52, q, x 






A, u(dyt) + Ajou;,(b,t) = AV3 A, ,U_; (t) Te2 






Ao u(a,t) + Aoou,(a,t) re A953 + AnyU_, (t) L1.3 






Figure Ll.1 - THE LPCM 





represents an approximation to an equation which has been linearized 
about an initial steady state, this portion should be zero if the 
model is to be valid, but cannot be assumed zero for any general 


model, Thus, u,(x) must satisfy equation L1.5 with boundary conditions 


given by L1.6 and 11.7. 


state u.(x). 


gorau] + = [P,(x)u, (x)] = 0 L1.5 


Ayy¥y(b) + Ay ots (b) = A,s L1.6 


Anzu, (a) + Agou, (a) aa A53 L1.7 


The second distribution implied by the model is the final steady 


the column after the step inputs and after all transients have had 
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This represents the lighter component concentration in 





time to die out. The final steady state u,(x) must also satisfy 
equation L1.5, but with the boundary conditions given by L1.8 and 


L1.9-. 
Ap u,(a) “- Angus (a) = Ao3 + Aon L1.9 


The third concentration distribution is the transient response 
u(x,t) - (Notes: the same notation is used for the transient 
response as was used in Ll.1 for simplicity) - which must satisfy 
equation L1.1, but with the boundary conditions L1.10 and Ll.11 
and the initial condition L1.12. In the interest of having a com- 
pletely defined problem and also from the practical constraint of 


proper model behavior, the final condition L1.13 must also be 


included, 
Ayju(b,t) + Apou,(b,t) = 0 eee 
Ao u(a,t) - Ant, (at) = 0 sens 8 
u(x,0) = u(x) = u, (x) - u,(x) Ligue 
u(x,co) = 0 L1.13 


The net result of considering these three spatial distributions 
is that the LPCM solution now is seen to be the solution of three 
separate problems: two ordinary differential problems L1.5 with 
different boundary conditions and one partial differential problem 


L1.1 with conditions L1.10 - L1.13. 
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The general usefulness of the LPCM must lie in the capability 
of using step response solutions to approximate the response to any 
arbitrary input functions of time. In utilizing the model in this 
general manner, the initial steady state u, (x)cannot be assumed to 
be zero and the model must be solved as the three separate problems 
described previously. 

For the purposes of presenting and solving the LPCM in terms of 
single step inputs, however, the assumption will be made in this 
thesis that Al3 and Ao3 are both zero and, thus, that the initial 
steady state distribution u, (x) is zero throughout the column, i.e. 
that the initial deviation from steady state is zero, Solution of 
the general problem implies solving for u, (x) in L1.5 and using it, 
along with u(x), in the transient boundary condition L1.12. 

The assumption that u, (x) = 0 for the purpose of evaluating 
the model reduces the solution to one partial differential problem 
Ll.1 in the transient response u(x,t) and one ordinary differential 
problem L1.5 in the steady state response u, (x) with boundary con- 
ditions given by L1.i4 and Li.15. Once the transient and steady state 


solutions have been determined, the total model solution u(x» t) for 


step inputs can be expressed as L1.16. Then u, (xs t) will satisfy 


the LPCM with A 3 = Ay = 0, The development of these solutions 


3 
begins with the application of the technique of separation of 


variables to the transient partial differential equation Ll.l. 
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u(x,t) = u,(x) + u(x,t) L1.16 


L1.2 SEPARATION OF VARIABLES APPLIED TO THE LPCM 


The method of separation of variables rests on the assumption 
that the solution u(x,t) can be separated into a product of two 


functions X(x) and T(t) as in eq. L1.17. This assumption is the 
u(x,t) = X(x) + T(t) L1.17 


Single most universally used analytical technique for solving 
partial differential equations. If this expression for u(x,t) is 
then substituted into the LPCM, the resulting expression is 


equation L1,18, where Ps and P, are defined in Ll.32 and Ll. 33. 
ee F - Ba 
PLXT + P3kT + PyAT = Xt L1.18 


If this expression is divided by XT,then the resulting expression 
in equation L1.19 represents a function of x equated to a function 
of t for all values of x and t. Thus, the two expressions must 
equal a constant - K®, 

The applicability of this technique to a given partial 


differential equation rests on the capability of separating the 


Bie tebe PRAT a L1.19 
> T 


resulting equation, as in equation L1.18. This step is, in general, 
not possible when the given partial differential equation is non- 
linear as are, for example, the continuous-spatial equation and 


polynomial-class partial differential equations (Davis, D-5, p.213). 
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The general condition for separability of a partial differential 
operator of the form of L1.20 has been shown by Murray-Lasso (M-5) 
to be the commutation of the operators Hy and Hi expressed in eq. 
Ll.2l. In this case Hes Hi, and Let are partial differential 


operators in x, t, and both x and t, respectively, In the case 
Ly ¢ (u(x,t) ] = Hy [u(x,t) ] + H, [u(x,t)] 11.20 
Hy [H, (u) ] = H, [H(u) ] L1.21 


of the continuous-spatial equation H. and Hy are defined by 


equations L1.22 and L1.23, where f(u,u,, uy) is a nonlinear function. 


H, [uj = gu 


Ll.22 


Hu] = f(uyu,u_) L1.23 


That these two operators do not commute is shown in eq. L1.24 and 
therefore the continuous-spatial equation is not separable, The 
~ [H,(u)] = le 35 - H, [a] — 

LPCM, however, does satisfy eq. L1.21 and is separable, as has 
been shown in eq. L1.19. 

The application of the separability assumption then results 
in the transformation of the partial differential equation in two . 
variables to a set of two ordinary differential equations, each 
of which depends upon the separation constant -K*, The two 
resulting ordinary differential equations are presented as the 


spatial equation L1.25 and the time equation L1.26. 
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T = -K9T L1.26 


L1l.3 SQLUTIONS OF THE TIME EQUATION 
The time equation L1.26 is an eigenvalue problem whose eigen- 


function solutions are easily found. Since u(x,t) in equation Ll.1l 
represents the transient portion of the desired solution, the 
separation constant has been chosen as -K® in order that each 
eigenfunction of the transient response approach zero for large 
values of time as in equation L1.13. The eigenfunctions which satisfy 
equation L1.26 are shown as functions of the eigenvalues “K,* in 


equation L1.27- 
ec) = exp(-K,7t) leq 


Ll.4 SQLUTIONS OF THE SPATIAL EQUATION 
The eigenvalue problem presented by the spatial equation L1.25 


and its associated boundary coriditions is much more difficult to 
solve than the time equation problem. In fact, statements about 
the existence, uniqueness, and continuity of solutions for this 
problem cannot be made in general (B-8). Statements of this nature 
require restrictions on the equation or the boundary conditions. 
There are some existence and uniqueness theorems for restricted 
cases presented in the literature (B-9) (K-1) (H-11), however, the 
existence, continuity, and completeness of eugentuniction solutions 
to regular Sturm-Liouville problems will be assumed and utilized in 
this section (B14) (I-2). 
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Several very Simple cases of the spatial equation L1.25 can 
be solved analytically. These simple cases are for Pi» Po and Py 
set equal to constants, and examples of these are solved in Chapter L4. 
In all of these simple cases and, in fact, for the problem L1.25 in 
general, the final steady state distribution in the transient solution 
boundary condition L1.12 is expanded in terms of the eigenfunctions of 
problem L1.25. For reasonably well behaved systems, such as a binary 
distillation column, it is expected that only the first few, say ten 
(10) or less, of the eigenvalues are really crucial to the transient 
behavior of the distillation column. The reason that only the 
smallest eigenvalues -K,* are important to the total transient 
behavior is that the time-eigenfunctions in equation L1.27 approach 
zero for very small values of time when the eigenvalues “KAS become 
larger. As mentioned in Chapter L4, the analytical solution eval- 
uations of the simple cases of the LPCM always required less than 
ten (10) eigenfunctions for three-decimal-place accuracy, 

These considerations offer promise of greatly simplified 
solution expressions for the LPCM if an analytical technique can be 
developed which offers a solution in terms of the first few eigen- 
values and eigenfunctions, AS mentioned in Chapter L5, there are 
numerous analytical and numerical procedures which offer this rep= 
resentation; however, conversion of L1.25 to an integral equation and 
then solution of the integral equation either by separation of the 
kernel (Appendix A2) or by approximating the integral equation by a 
system of linear algebraic equations (Appendix A3) will be investigated 


in this thesis, 


Seo 





The integral equation technique consists of converting the 
differential equation to the Liouville Normal-Form equation and then 
converting the Liouville Normal-Form equation to a Fredholm II- 
integral equation. The first conversion procedure begins with a 
test for self-adjointness and a transformation of L1.25 to a regular 
Sturm-Liouville equation L1.40. The end result of these manipula- 
tions is the set of spatial equation eigenvalues K, in equation L1.49 
and the set of eigenfunctions W(z) in either equation A2.15 or 
A3.10 with the transformation relations L144, L1.45, L1.37, and 
L1.39 leading to the eigenfunction solutions X (x). The following 
sections present the development of these transformations and the 
equation solution procedures, 

L1.5 SELF-ADJOINT CONDITIONS AND TRANSFORMATIONS 

If the spatial equation L1.25 can be shown to be self-adjoint, 
then the spatial eigenvalue problem is a regular Sturm-Liouville 
eigenvalue problem. If the spatial equation is not self-adjoint, 
then the variables can be transformed into a self-adjoint equation 
which then represents a regular Sturm-Liouville problem. In an 
analytical analysis such as employed in this chapter, it is essential 
that at some step the equations be expressed in a‘format for which 
real eigenvalues and complete eigenfunctions are guaranteed, The 
regular Sturm-Liouville problem is such a format (B-14) (I-2). 

The necessary and sufficient condition that the spatial equation 


L1.28 be self adjoint is expressed by equation Ll.29, When this 
€® § ~~ / 2 = : 
P (x)X “ Po(x)X “ [P, (x) + K?] X= 0 L1.28 


P, (x) = P,(x) L1.29 
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condition is applied to the coefficient polynomials in equations 
Ll. 30, L1.31, L1.32, and L1.33, the polynomial coefficient relations 
L1.34 and 11.35 result. The self-adjoint form of equation L1,28 


is then given by equation Li. 36. 


ntl al 
P - eee A 
(x) = Psa Pa L1. 30 
n+l jel 2 5 
Po (x) 1 15-1” aT hale: 
Pa(x) = q +5 (25p. + .x)xo7t Ll. 32 
3 Io yt 35 J : 
P,, (x) +P (sn) Typ. + Jol 
yl) = a) + Ft LIP. * jy Xe 233 
Q523 + JP, =) J iae 15 92 <50n Ll. 34 
an = 0 L1. 35 
d_ PP, (x)X] + P,(x)X = -K2x Ll. 36 
re fa) 4 


If the spatial equation L1.28 turns out to be non-self-adjoint, 
then the transformation relations L1.37, L1.38, and L1.39 (M-2) 


result in the self-adjoint equation L1.40. This transformed 


P(x) = exp [jf P3(x) ax] 11.37 
xP LS, = pe) 
Q(x) = Py(x) , P(x) L1. 38 
P, (x) 
R(x) = P(x Ll. 39 
ate) 
a. [P(x)X} + Q(x)X = -K®R(x)x L1.40 


equation L1.40 is then seen to be of the same form (Sturn- 


Liouville) as the self-adjoint equation L1.36. The spatial 
~116- 





eigenvalue problem thus reduces to a regular Sturm-Liouville problem 
with the differential equation expressed by L1.40 and the boundary 


conditions by L1.41, L1.42, and 11.43. 


A, X(b) “: A, 2X(b) = 0 L1.41 
An, X(a) + A, X(a) = 0 L1.42 
Aj Ago - Ay aAno #0 L1.43 


The next step in the solution procedure is to convert the 
regular Sturm-Liouville equation L1.40 to the Liouville Normal-Form 
equation L1.46, Equation L1.25 could have been converted directly 
to an integral equation but this was not done for several reasons. 
Firstly, the analytical manipulations in the conversion procedure 
did not lead themselves to the approximation techniques so necessary 
for finding the final solution. Secondly, and much more importantly, 
there is no guarantee that a direct transformation would produce an 
integral equation with a continuous kernel. In fact, in most cases 
it would not. If either the kernel separation method in Appendix A2 
or the linear algebraic method of Appendix A3 were employed for 
solving an integral equation with a discontinuous kernel, there would 
be no guarantee that the eigenvalues would be real and distinct. Such 
& guarantee can be given if the kernel is continuous and symmetric. 
L1.6 TRANSFORMATION TO THE LIOUVILLE NORMAL-FORM EQUATION 

Transformation on both the dependent and independent variables 
of equation L1.40 can be used to greatly simplify the regular Stura- 


Liouville equation (B-14) (I-2). If new variables W and z are defined 


be 





as in equations L1.44 and L1.45, then the equation L1.40 is transformed 


to that of equations L1.46 and L1.47, where the primes are derivatives 


X = W L1. 44 
4 PCx)RCX) 
Z = Wy, R(x) dx L1.45 
P(x 


with respect to x and the equations are expressed as functions of z. 


aw 
dz* + [d = M(z) JW = 0 L1.46 


1 lal 1 2 
(2) a +) +2@) 
=~ L*2 
L (5 F(R) - LG) }+Q 11.47 


Whereas the original equation L1.40 was defined over the interval 


a. 


1A 


x <b, the transformed equation L1.46 is valid over the interval 


0 


1A 


zZ<c, where c is defined in equation L1.48 and A in equation L1.49. 
b 
c= f Y BGY oe L1.48 
a 
P(x 
A = K? L1.49 


The boundary conditions of L1.41 and Li.42 are also transformed to 
the new boundary conditions on W(z) expressed in equation L1.50 


through L1.56, where the dot is d/dz and the prime is d/dx. 
| : * ° 50 

D,,W(c) + D, ,W(e) 0 LaenS 

Dz, W(0) + Dop¥(0) = 0 Ll. 51 


1 i 
Dyy = Ayy mee + me) 
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P(b 
a 


it 
Do, ™ Ap, - Aae es res 32) L1. 54 
ay «\Pta R(a 
Doo = Ano R(a) B55 
P(a, 
Dy 4255 - Dy oDo1 #0 Ll. 56 


L1.7 SOLUTION OF THE LIOUVILLE NORMAL-FORM PROBLEM 

Now that the spatial differential equation L1.28 has been 
transformed to the Liouville Normal-Form, the next step in the 
solution of the LPCM consists of finding the eigenvalues A; and the 
eigenfunctions W, (2) for equation L1.46 with the boundary conditions 
L1.50 and L1.51, As discussed in subsection L1.4 the technique to 
be used is that of converting L1.46 to an integral equation and then 
solving the integral equation, 

The conversion procedure is presented in Appendix Al and two 
alternative solution procedures are presented in Appendices A2 and 
A3. The separable kernel procedure presented in Appendix A2 could 
be employed for specialized cases of the LPCM. The solution procedure 
in that case would be to approximate the kernel Al.17 by a finite 
series A2,2 and then to apply the techniques of Appendix A2, This 
would involve a great deal more analytical analysis for any given 
LPCM, but it is anticipated that the resulting solutions would be 
more accurate than those obtained by the methods of Appendix A3, 

The relative merits of these two solution procedures in terms of 
solution accuracy and computation time would have to be examined for 


any given LPCM, 
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The second integral equation solution procedure in Appendix A3 
consists, in essence, of approximating the integral by a finite sun, 
The procedure is easily programmed on a digital computer and is a 
very direct method. For these reasons, the procedure in Appendix A3 
will be used to evaluate the LPCM in this thesis, 

Once the eigenvalues K, and eigenfunctions X(x) of the spatial 
equation have been found, the next step in finding the total solution 
to the LPCM is to ensure that the initial condition L1.12 is met, 
This requires that the initial distribution u(x) be exproeees in 
terms of those eigenfunctions X (x). 

L1.8 THE EIGENFUNCTION BXPANSION 

The eigenfunction expansion begins by expressing u, (x) in equation 

11.12 as a sum of constants c, multiplied by each eigenfunction X, (x) 


as in equation L1.57, when n is the number of eigenfunctions, Since 
> 1 
(2%) $_, 04%4 L1.57 


the integral equation solution technique only gives n eigenfunctions, 
then equation L1.57 can only be satisfied for n points x. in the 
interval a< x <b. Thus, equation 11.57 must be discretized to 


equation L1,58. Equation L1.58 then represents n-simultaneous equations 
(x.) = £_0,X;(x,) L1. 58 
U(x; £_°4 j (Xs : 


in neunknowns which can be expressed as the matrix equation 11.59. 


The solution to L1.59 is then given by L1.60. Once the matrix C has 


U, = XC igo? 


Ge XU L1.60 
e =e 





been found, the total solution can be expressed in terms of the 
transient and steady state portions. 
L1.9 THE TOTAL SOLUTION TO THE LPCM 

The total solution to the LPCM can now be expressed in terms of 
known functions u(x) and u(x,t). The steady state problem L1.5 
with boundary conditions L1.14 and L1.15 is a standard linear boundary 
value problem for which solutions and, in fact, computer subroutines 
exist. The steady state solution is determined in this thesis by 
Scientific Subroutine Program - LBVP in reference (I-3) and is 
described in greater detail in Chapter L2. 

The major portion of the LPCM, for which computer programs have 
not been developed, was the solution to the transient boundary value 
problem. For this reason, this portion of the LPCM received the 
greatest portion of solution effort in this thesis, It is expected 
that the steady state equation couid be soived by the same integral 
equation techniques (leading to a final matrix inversion rather than 
an eigenvalue problem) as were used for the spatial equation portion 
of the transient response. The conversion would be more involved 
because the boundary condition constants Aaa and Aon in equations 
L1.14 and L1.15 are not zero. The conversion of both the steady state 
equation and the spatial equation of the transient response, and the 
incorporation of the solution technique for both of these problems 
in the same program would undoubtedly result in a great savings of 
computation time over that spent using the LBVP subroutine and the 
integral. equation techniques separately. This could only be determined 


from further studies. 
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The efforts of the previous sections of this chapter have led 
to the transient solution to the LPCM due to step inputs. The total 
solution was given in equation L1.16; the transient portion of that 


solution can now be given as equation L1.61. The total solution to 
n 
= - “ 
u(x,t) = e caky (x) exp(-K, *t) L1.61 


the LPCM is then given by equation L1.62. 


u(x,t) = u(x) + Z 1 04X,(x) exp(-K, ®t) ER 





The next chapter (L2) develops a digital computer program for 
performing these manipulations and calculating the solution to a 


special version of the LPCM, 
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CHAPTER L2 


THE PARTIAL LINEAR=POLYNOMIAL BOUNDARY VALUE PROGRAM (PLPBV) 


The very general solution technique developed for the LPCM in 
Chapter Ll is based upon mathematical transformations whose analytical 
expressions, for general nth order polynomial coefficients, are too 
complicated to use in a numerical solution evaluation. This chapter 
presents an application of the integral equation technique of 
Chapter Li to the first-degree polynomial model (PLPBV) and a digital 
computer program for evaluation of the solution. 

L2.1 PRESENTATION OF THE MODEL (PLPBV) 

The Partial Linear~Polynomial Boundary Value (PLPBV) model is 
presented in Figure L4.1. ‘This model is a subcase of the LPCM of 
Figure L1.1 in which P,(x) is a linear polynomial and P, (x) isa 
constant. This model represents one step upward from the simplified 
LPCM examples of Chapter L4 and, at the same time, the simplest form 


of the LPCM which still has a spatial varying coefficient. 





2 [P,(x)u] + 9_ [Po(x)u] = du 
ie ie 





Where: wu = u(x,t) 





P,(x) = By PB, > 0 





Po(x) = a + Q xX 





Auld, t) + Ayou,(b,t) = Ay y0_, (+t) 






Aula, t) + Anou,(a,t) > A5),U_, (+) 






Figure L2.1 = THE PLPBY MODEL 


EIDaS: 





The development of the solution to this model now parallels that of 
Chapter Li for the general LPCM except that the transformation 
relations can be expressed specifically in terms of the polynomial 
coefficients. 
L2.2 THE TRANSIENT AND STEADY-STATE EQUATION FROM PLPBV 

The first step in solving the PLPBV model is the expression of 
the separate transient and steady-state problems as in subsection Ll.1l 
and equation L1.16. The transient and steady-state problems are 
presented in detail in Figure L2.2, where u,(x,t) is the total solution 


to PLPBV. 
u(x,t) = u(x) + u(x,t) 


Transient Portion: 


PU + Pau. “ Pu = UL 


where 3 PS 2D. 

P. 1. + Q,x 

Ei SS) 
A, uld,t) + A, ou,(b,t) = 0 
A,.u(a,t) - Annu, (a,t) = 0 
u(x,0) = <u_(x) 
u(x,co) = 0 


Steady-State Portion: 


Pius “ Pou. “ Pu. = 0 


where; PS = Po 


P3 = oi + q4 x 


Pye 9) 





(Continued on next | 
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Ajyug(b) + Ayous(b) = Ayy 


A,.u(a) “ A,ot(a) = Aa, 


Figure L2,2 - TRANSIENT AND STEADY-STATE 
PROBLEMS OF PLPBV 





L2.3 SOLUTION OF THE STEADY STATE PROBLEM 

The solution of the general LPCM steady state has been discussed 
in subsection L1.9. The subroutine LBVP, Linear Boundary-Value 
Problem, and the associated subroutines applicable directly to the 
general LPCM are presented in Appendix A5 for completeness and for 
the convenience of the reader, Subroutines LBVP and GELG are taken, 
less comment cards, directly from the literature (I-3), whereas 
subroutines AFCT, DFCT, and FCT have been written specifically for 
the general LPCM. The user must supply his own output subroutine 
OUTP as per reference (I-3). 

Once LBVP is called in the main portion of PLPBV tnen it calculates 
the steady-state values us (Xs) for use in equations L1,58, for the 
final transient eigenfunction expansion, and in L1.62, for the total 
solution evaluation at the discrete points X se The user OUTP sub- 
routine must be set up to return u(x.) values to PLPBV at exactly 
the desired X 5° 

Some comments concerning the efficiency of LBVP in the applica- 
tion must be made, It is obvious that since both FCT and DFCT return 
zero values that a simpler version of LBVP for this special case 


could be writtenar possibly found in the existing literature by 
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further searching. As mentioned in L1.9, the integral equation 
techniques might also be applied here, 

For the initial purposes of evaluating the LPCM and testing the 
transient portions of PLPBV, the steady state constants were fed in 
as Column 9, Matrix 3, in Appendix A4Y, Subroutines LBVP, AFCT, DFCT, 
and FCT have all been compiled, tested, and shown to operate properly. 
L2.4 SOLUTION OF THE TRANSIENT PROBLEM; ANALYTICAL ANALYSIS 

The solution to the transient portion of PLPBV defined in 
Figure L2.2 will be obtained by applying the transformation rela-= 
tionships to L2.5 to convert the spatial equation equivalent to L1.25 
to the Fredholm II-Integral Equation Al.21. Tnis analysis will be 
divided into an anaiytical anaiysis in which the specifics of the 
transformations will be presented and a computational analysis in 
which subroutine PLPBV will be explained and presented. 

The analytical analysis of the transient portion of PLPBV begins 
with separation of variables, The time equation eigenfunctions 
resulting from the application of the separation-of-variables 
technique are the same in PLPBV as they were for the general LPCM 
in equation L1.27. The resulting spatial equation is given below as 


equation ee Te 


me + Pk Tae kK | Kd L2.7 
where s Ps a 
Ps ae “ q, x 
P= 
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The next step in the analysis of the spatial differential 
equation L2.7 is a test to see if the equation is self adjoint or not, 
If it is self adjoint, then no transformations are necessary because 
then the equation can be written directly in the form of Ll. 36. 
Application of the self-adjoint test conditions Ll.34 and L1.35 
results in equation L2.8 which shows that equation L2.7 is not self 


adjoint and that transformations wiil be required. It is interesting 
q,=0 , but q, #0 in PLPBV L2.8 


to note at this point that aq, = O is satisfied for the simplified 
version of the LPCM solved analytically in Chapter L4 showing that 
it is already self adjoint. 

Since equation L2.7 is non self adjoint, the transformations 
Ll. 37, L1.38, and L1.39 can be applied, as in equations L2.9, L2.10, 
and L2.11, resulting in the self-adjoint equation L1.40. With these 
functions defined, the second transformation, to the Liouville Normal- 
Form equation, can be utilized as in equations L2.12 and L2.13 


resulting from equations L1.44% and L1.45. Application of transformation 


P(x) = exp [(x~a) a,/p, + (x?-a)q,/2p,] 12.9 
Q(x) = P(x)a,/p, L2.10 
R(x) = P(x)/p, ea 
X(x) = W(x) {BS exp [-(x-a)q,/2p, - 
(x3~a?)q,/4p, | L2.12 
2= (xea)/ JP, bei 3 


12 7= 





relations L1.47, L1.48, and L1.49 results, after some manipulations, 
in equations L2.14, L2.15, and L2.16, Inverting L2.13 and substitut- 


ing for x in equation L2.14 gives M(z) in L2.17. 


M(x) = do* + 2a4)x + 4,7x? 3a, 


L2.14 

4D, “a 
Cc = (b-a)/P, L2.15 
X — K? . 2 ako 


M(z) = 4,° + 24,0, ( vB5 2 + a) + 4,7 ( vB 2 + a)? 


“ 39, 
3 L2.17 
The next step in the trensformation to the Liouville Normal- 
Form equation is the evaluation of the boundary conditions in L1.50 
through L1.55. These relationships are given in equations L2.18 


through Ze 22 e 


t= Ay ola. ¥ q, 0) /2p, L2.18 
Dia = Apa/ V% L2.19 
Day = Ap, 7 Analy + 938)/2P, | eee 
Dog = Ana! V Py L2,21 


DyDog = DypPoy = (Ay1490 ~ ApyAro)/ VP, 


P | 
~ Ay pAondy (b-a)/2p,? # 0 E2522 
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The final transformation which must be made is the conversion 
of the Liouville Normal-Form equation to the Fredholm II-Integral 
Equation following the steps in Appendix Al. However, all of the 
functions and constants in equations Al.20 through Al.24 are now 
known in terms of PLPBV, so the transformation is complete. With 
these known transformation relationships and with the solution 
technique of Appendix A3, a digital computer subroutine can be 
written for computation of the solutions Xx) the eigenfunction 
expansion, and the computation and plotting of the final solution 
u (x,t). 
L2.5 SOLUTION OF THE TRANSIENT PROBLEM; COMPUTATIONAL ANALYSIS 

Once a complete analytical analysis of any specific subcase of 
the LPCM has been accomplished, the next step is to write a computer 
program to evaluate the solutions utilizing analytical transformations 
and integral equation solution techniques, Such a computer program 
will require a significant amount of work to write and computation time 
to check out. However, it is expected that, in the end, the program 
Will require far less computation time than currently available programs 
which solve the discrete~plate equations, It is expected that a 
further significant savings in time would result from the design of a 
special purpose computer to solve a model of sufficient degree to 
attain desired accuracy for a specific column or type of column. 

All computer programs utilizing the methods of this thesis will, 
by the nature of the method, follow a certain format. This subsection 


presents a description of the computation steps, in flowchart forn, 
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which this solution method uses, This flowchart is presented in 
Figure L2.3, and a Fortran IV program listing the details of the block 
steps in Figure 12.3 for the special case of PLPBV is presented in 
Appendices A5 and A6, 

The PLPBV subroutine presented in Appendix A6 is intended as a 
description of the steps necessary for solving LPCM problems. Sub- 
routine PLPBV was compiled and executed for several very simplified 
cases but would require a great deal of work to be useful in general, 
It is estimated that PLPBV would require between 100 and 200 man-hours 
of programming time and 1 to 2 hours of computation time to perfect, 
to consider all special cases, and to debug, The primary benefit 
from this labor would be the fact that PLPBV only requires about 15 
seconds of computation time, and thus, complete column solutions using 
PLPBV could be generated in less than one minute for 20 spatial points 
and 9 times. 

This chapter has presented a detailed description of the steps 
necessary to practically apply the general integral equation solution 
technique developed in Chapter Li to a special subcase, PLPBV, of the 
LPCM. These steps are described analytically by applying the trans- 
formation relations to PLPBV and are described numerically by a flow- 
chart of the necessary computation steps and by a very basic Fortran IV 
digital computer subroutine showing some of the details of the computa- 
tion steps, The next chapter describes the steps necessary to determine, 


the LPCM for a specific distillation column, 


moO = 
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CHAPTER L3 


DETERMINATION OF THE LPCM FOR EXAMPLE DISTILLATION COLUMNS 

The Linear Polynomial-Coefficient Model (LPCM) has been developed 
and solved in this thesis as a continuous-spatial model of the transient 
composition behavior of a binary plate distillation column. The purpose 
of this chapter is to show the steps required to determine the LPCM 
for given distillation columns and to discuss the approximations 
involved in developing the LPCM, Two columns having the same feed and 
output compositions but different reflux rates and numbers of plates 
are used as examples to show the development of the LPCM. This chapter 
represents an equating of Figure M4,2 to Figure Ll.1 in terms of the 
two example columns. 
L3.1 GENERAL STATEMENT AND DISCUSSION OF THE MODELING STEPS 

The first step in the determination of the LPCM for any distilla- 
tion column is to determine the slope of the equilibrium curve, m(x), 
as a continuous function of x. The function m(x) is defined in equation 
M4.5 and is presented here in detail as equation L3,1. Based upon the 

m(x) = df(u) ~ ed eS pal 
Qu [1 + Ca-1)u, (x)? 
u=u, (x) 

initial steady-state us (x5) represented by the McCabe-Thiele diagram, 
the function m(x) can be expressed exactly at only a discrete number 
of points. The major approximation involved in determining the LPeM, 
aside from the original linearization of the CS, is the expressing 


of m(x) as an n-th degree polynomial in x, as in equation L3.2. 
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n 
m(x) = mx =m, + mx + Mx* + ues LBFZ 


The continuous-spatial steady-state model has not been explicitly 
emphasized in this thesis because the direct approximation of the 
oo steady-state equilibrium curve slope seems to offer a more 
direct approach. The continuous steady-state model is equivalent to 
Figure M2.1 with all time derivatives set equal to zero, If this model 
were solved for the continuous steady-state u; (x), then the continuous 
m(x) could be found from equation L3.1. This is impractical for the 
LPCM determination because the approximations necessary to develop and 
solve the continuous steady-state model introduce more error than the 
direct determination of m(x) from the discrete steady-state model. 
Thus, continuous steady-state models are not used in this thesis and 
are very seldom found in the literature, 

The accuracy of the LPCM is expected to depend directly upon the 
accuracy of m(x). Theoretically m(x) can be determined to any desired 
accuracy uSing a polynomial in equation L3.2 of sufficiently high 
degree and using either point-by=point or leastesquares techniques. 
The integral equation solution technique of Chapter Ll is valid for 
any degree polynomial, but as the degree gets higher than one the 
transformations become overwhelmingly complicated and the programming 
time necessary to implement the transformations and consider all of 
the special bases becomes an order of magnitude larger, Surprisingly 
enough, the computation time, once the proper programs are written, 
is not expected to increase significantly over that for PLPBV in 


Chapter L2 because the basic steps of Figure L2.3 remain nearly the 
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same and all of the special cases are parallel paths for computation. 
It is expected that models of degree one or two will prove to be 

, adequate to specify the column behavior sufficiently for the purpose 
of control. 

Once the continuous function m(x) has been determined, the model 
derivation is essentially complete. All that remains is to calculate 
the polynomial coefficients Ps and qs and the boundary condition 
constants ae The complete expressions for evaluation of these 
constants in terms of a second degree m(x) are given in Figure L3.3. 
L3.2 APPLICATION OF THE MODELING STEPS TO TWO EXAMPLE COLUMNS 

The McCabe-Thiele diagram for a five-plate binary distillation 
column is presented in Figure L3.1 and the corresponding diagram for 
an eleven=plate column is presented in Figure 12.2 in Section 1(I). 
These two diagrams have been designed to have the same input and output 
compositions to simplify the calculations and to provide for easy 
comparison. The discrete values of m(x) have been calculated using 


equation L3.1 and are presented numerically in Table L3.1 and Table L3.2 


and presented graphically in Figure L3.2. 


5 OFe2 0.94 


Table L3,1 = THE DISCRETE m(x) FOR A FIVE PLATE COLUMN 
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Figure L3,2 - EQUILIBRIUM CURVE SLOPES, m(x), FOR EXAMPLE COLUMNS 
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Model Equation 
g?_ LP; (x)u | +O_ LP, (x)u ] =du 
Qx* ax gt 


Where; 


P, (x) =P, + BX + Dox? = [(B +m) + mx + myx? ]/2 


Thus: p, = (B+ m)/2 
p, = m,/2 
i M 
Po = m/2 
P5 (x) = 4) + a,x + q5x* = [ (B - n.) ~ mx - mx* | 
Thus: als : (B - nm.) 
a 
Qo al in) 


Both P, (x) and P, (x) must be expressed for both the upper 

and lower sections of the column using BY and B, - 

Boundary Conditions 

Upper Ay, = - m(x) + 7 | 
"| 
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b=1.0 dx 
a=X_p x= 1,0 
Ayy = Mo 7 mh + 7 
me 6 
sath 
Any - “[3,, es ~ ae + JVI 
fea, ao LE 2% £ 
eo 
Aoi = -F/V 
Lower ~ eg! + 2moXn) 15 oY 
beXxX- us =iji- ~ M,Xp ~ My 
: . : ~ Mo 
f- 2 ot 
Ae3 8 5 
Aol 


Feed Condition F/V = (B, - B )/4 


Pigure L3.3 - THE COMPLETH LPCM WITH SECOND DEGREE POLYNOMIALS 
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Table L3.2 = THE DISCRETS m(x) FOR AN ELEVEN PLATE COLUMN 





The next step in the modeling procedure to develop the LPCM for 
these two columns is to determine the coefficients of a polynomial, 
m(x), of sufficient degree to meet accuracy requirements based upon 
the discrete points given for m(x), In the case of the 5-plate column 
a polynomial of degree 5 could be found which would exactly match each 
of the 6 points given, and for the ll=-plate column a similar polynomial 
of degree 10 could be found. However, for both of these columns a 
2nd-degree polynomial is expected to be accurate enough. ‘Two second 


degree polynomials have been found as examples using the points, 


i 


x = (0.0, 0.2, 1.0), for the 1l-plate column and the points, 


if 


x = (0.0, 0.6, 1.0), for the 5-~plate colum. The resulting m(x) 
polynomials are given by equation L3.3 for the 5-plate column and 


Oe 





L3.4 for the 1l-plate column. 


m(x) = 2.31 - 4.17x + 2,22x? 5-plate oes 
m(x) = 2.31 - 5.32x + 3.37x? li-plate L3.4 


Once the function m(x) has been determined to the desired accuracy, 
the complete LPCM for the column can be calculated from the equations 
of Figure 13.3. The 2nd-degree m(x) polynomials in equations L3.3 and 
L3.4 result in the complete LPCM for the 5=plate column in Figure L3,4 
and the complete LPCM for the li-plate column in Figure L3.5. It should 
be emphasized that the approximations and calculations used to obtain 
equations L3.3 and L3.4 were chosen to greatly simplify the calcula- 
tions. Normally, the polynomial coefficients for m(x) should be deter- 
mined specifically for the region of application of the model equation 
instead of using an overall m(x) as was done in equations L3.3 and 
L3.4. The column may also be modeled using more sections than the 
upper and lower sections used in this thesis. 

L3.3 .GENERAL COMMENTS AND SUMMARY 

The real utility of the LPCM results from the fact that as the 
number of plates in the column increases, the model complexity remains 
nearly constant for realistic approximations of m(x). The computation 
time and complexity of discrete models, however, increases rapidly as 
the number of plates increases because the dimensions of all of the 
matrices in the model increase with the number of plates, The LPCM is 
valid for components with non-constant relative volatilities as well 
as for constant relative volatilities, but only the comstant relative 


volatility equilibrium curve has been used in this thesis, Also, the 
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5-PLATE COLUMN 
Model Equation 
92 [y(x)u] +0 [P,(x)u] = du 
gx? gx ot 
Where: Bis apt Py x oF Pox® 
Po x) = Qa, + a,x + Qox* 
Upper Section Lower Section 


Py = -2-09 a) = 4.17 Py = 72.09 gq, = 4.17 


Boundary Conditions 
A, u(d,t) + A ou, (d4*) =A t Ay U_, (t) 
An, u(a,t) + Ajou, (ast) = Ang + Apy U_, (t) 
Upper Section b=1.0, a=0.5 
Aa = 0,91 Ajo = 0, 36 Ay 3 = 0,0 Aaa = (0506 
Ady = ys Ano SO Pe Ag = 0,0 Aon = -0,17 
Lower Section b=C. 5, a=0,0 


21 
Feed Condition 


22 


F/V = 0.17 


e L3.4 = THE COMPLETE LPCM FOR THE 5=PLATE COLUMN 





LPCM can be used for multiple feed and sidestreams and can be applied 
in any manner desired to smaller sections of the column, Each sube- 


section of the column can then be modeled using the LPCM, 
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l1-PLATE COLUMN 


Model Equation 


2 [P,(x)u] + 0_ [P,(x)u] = Qu 
ga Cad +O a 





Where: Bayt P, x + Pox 


nfs 
DeQmetde + ox + qox? 


Upper Section 0.5=5xS$1.0 Lower Section 0.0 Sed, D 


= 1,40 q 1.92 Cn -0.79 
Beer oo) 66 a -2.66 4) = 5.32 
= 1.69 qo =3.35/ 1.69 Qo = =3, 37 


it owt 
Wn 
® 
WW 
nN 
ae) 
bt 
bo ou il 


Boundary Conditions 
A, u(b,t) + Ayou,(b,t) = Ay3 + Ay U_, (t) 
An, u(a,t) + A, ou hs t) = Ang + Aon U_, (t) 
Upper Section b=1.0, a=0.5 
Any = =9,07 Ago = 0.51 A3 = 0,0 Ann = «1,03 
Lower Section b=0.5, a=0.0 
mya = 2493 Ayg = OSL Ay = 000 Ay = -2.03 


peel, 51 4A, = 1.52 3 = 0.0 Aon = 0.00 


21 
Feed Condition 


oe 


R/V - i053 


- THE COMPLET“ LPCM FOR THe J1-PLATE COLUMN OF FIGURE 12,2 





The two models in Figures L3.4 and L3.5 are not solved in this 
thesis. Their solution will require the development of a computer 
program. based upon the integral equation solution technique of Chapter Ll 
and using the computation steps of Chapter L2. This will require an 


ari = 





extensive amount of analytical analysis, computer programming, and 
testing. The result of this labor is expected to be a general, very 
fast, solution program for the transient behavior of a large class 

of distillation columns which describes the column behavior accurately 
enough for control applications, Extensive research would be required 
to justify this, however. 

This chapter presents the steps necessary to derive the LPCM for 
general distillation columns, The details of these steps are then 
applied to a 5-plate and an ll-plate column and the resulting LPCM’s 
presented. The next cnapter solves analytically and computes numeri- 
cally the solutions to some simplified LPCM's in which m(x) = Moe a 


constant. 
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CHAPTER L4 


ANALYTICAL SOLUTIONS TO APPROXIMATED COLUMN EQUATIONS 
The purpose of this chapter is to analytically solve and numer- 


ically calculate the total solutions to four simple subcases of the 
Linear Polynomial-Coefficient Model (LPCM) defined in Figure Ll.l. 
The reason for doing this is to provide some examples of simpler 
solutions upon which to base an understanding of solutions to more 
sophisticated versions of the LPCM. The procedure to be followed 
will be to define a simplified version of the LPCM, state the 
problems to be solved, present the analytical solutions to the 
problems, present a computer program used to sum the series, and, 
finally, to present graphically the simplified model solutions. 
L4.1 A SIMPLIFIED LPCM 

A greatly simplified LPCM is presented in Figure L4.1 by 
equations L4,1 through L4.5. This model represents one step below 
the more complicated model presented in Chapter L2, This simplified 
LPCM appears in many places in the literature, only it is known by 
the different names and uses listed in Table L4.1, for a representative 
sample of the literature, This model is usually solved by analytical 
methods in the literature and used as an example of the application 
of the method of separation of variables. 
L4.2 APPROXIMATE EQUATIONS TO BE SOLVED 

There are four specific problems involving the simplified LPCM 
which will be solved analytically in this chapter. These problems 
are presented in Figure L4.2 and Table L4.2, Two of the problems 
involve the solution for the top step response of the heat equation 


altro = 





for the two cases of zero bottoms composition and zero bottoms outflow. 
The other two problems solved are the Taylor Diffusion Model (TDM) 
equation for positive and negative values of the center constant Qo° 

A visualization of what is happening in terms of distillation 
column concentration or, by analogy, temperature in a rod is presented 
in Figure L4.1. If the response desired is the total response to a 
step change in feed composition, then the equations solved in this 
chapter apply directly to the bottom half of the column in the range 
O<x<¢1, where x = 1 is the location of the feed tray. The heat 
analogy would be to visualize a rod with length in the range 0 <x<1l 


with a step change in temperature at the top. 


3 Uu +0 Uu Zidy 
goa [Pot a Lao! a 


P, (x) = Py 

P(x) = 4, 

Ayyu(b,t) + Ayouy(d,t) = AqyU_y (t) 
A, ,ula,t) + Aoou. (ast) 30 


Figure L4,1 - A SIMPLIFIED LPCM 
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AUTHOR(S REFERENCE EQUATION NAME OR USE 


Crank C8 Diffusion Equation 









Gould 





Taylor Diffusion Model (TDM) 





Jackson and Pigford J-l Transient-Diffusion Equation 





Lee 









Pollock, Brown » Convective Transport Equation 
and Dempsey Pe5 







Stone and Brian 






Bennett and Meyers 







Boyce and DiPrima 


Berg and McGregor 









Jury 
Heat Equation 






Powers 






cagan 






Tsang 


Woodle 







Table L4,1 - LITERATURE NOMENCLATURE AND USAGE 
GF THE SIMPLIFIED LPCM 
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Column or Rod Column Concentration Chart 


Transient Distributions 


Visual Graphical 


Figure L4.2 = VISUAL AND GRAPHICAL SOLUTION FORMAT 


Equation Boundary Conditions 
i Po % "11 “Mio "Azi “Abe) Aly 


Table L4,2 - LPCM PROBLEMS TO BE SOLVED ANALYTICALLY 
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L4.3 ANALYTICAL SOLUTIONS TO THE SIMPLIFIED LPCM EXAMPLES 

This subsection develops the irfinite series representations 
of the analytical solutions to problems UAL, UA2, UB1, and UCl which 
have been summarized in Table L4.2. The manipulations leading to 
these solutions represent no more than a series of mathematical 
exercises. The details of these solutions are included in this 
thesis for the sake of completeness and for the convenience of a 
reader who may want to relate the solution techniques used for the 
general LPCM solution in Chapter Ll to the much more familiar 
techniques used in this chapter. 
L4.3.1 SOLUTION OF UAL 

The total solution to UAl is given by equation L4.1; the 


development is presented in the steps below, 


Statement of Problem: Step L4.1 
x 
= u=U_; (t) 
a u t 0 1 
u(1,%) ~ U_, (+) 
u=0 
u(o,t) = 0 
0 = 
u(x,o) = 0 
Transient plus Steady State: Step L4.2 
u(x,t) = u(x) -: up(x,t) 
cha: = Q a7un, = Unt 
u. (0) «= 0 up(o,t) = 0 up(x,0) = <u .(x) 


u.(1) =1 up(1,t) eo Up Xse0) = 0 
Steady State Solution: Step L4,3 


u(x) = xX 





Transient Solutions Step L4,4 
up(x,t) = X(x) * T(t) 
T, (kK) = exp (-K?t) 
X (x) = Ayeos K x + A, sin K x 
Un(o, t) = 0 —- A= 0 
up(1,t) = 0 defines eigenvalues sin K = 0 
gives K = nnn . 
Kigenfunctions ares 
a, sin (nqx) exp (-K#t) 
Expanding to meet the initial conditions 
Up(x, 0) = -u_(x) = <x Tm, sin (n7x) 
Using the orthogonality property of eigenfunctions: 
rr (=x) sin (myx)dx = Za, fe sin (mpx) sin (nqqx)dx 


m 
(-1) = an 


My Vs 


EKigenfunction constants: 


n 
as 2(-1) 


ney 


Transient Solution: 


fos) 
up(x, t) = 2 z (71) sin(ngx) exp [-(nna)?t | 
= ner 
Total Solution to UAL: 


u(x,t) = x + 2 Z_, (-1)"sin(nex) exp [=(neree) 2t 7] L4.1 
Tey 
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L4.3.2 SOLUTION OF UA2 
The total solution to UA2 is given by equation L4.2, the 
development is presented below, 
Statement of Problem: Step L4.1 
ayy, 7 Uy = 0 
u(1,t) = U_,(t) 
u,(o,t) = 0 


u(x,0) = 0 





Transient plus Steady State: | Step L4.2 
u(x,t) = u(x) -- u(x,t) 


Bete O Cty = Uy 
u_(1) =] un(i,t) = 0 up(x,0) = -u, (x) 


Us,(0) = 0 Up, (0,t) = 0 unlx,2) = 0 
Steady State Solution Step L4.3 
u, (x) = 1] 
Transient Solution: Step L4.4 


Un(xst) ‘- X(x) - T(t) 

S(t) = exp (-K*+) 

X (x) = A, cos x x + A, sin <x 

Un, (04+) = 0 gives Ap = 0 

up(1, t) = 0 defines eigenvalues cos K = 0 
gives Kk. = nny a 

BEigenfunctions are: 

&, GOS (nx) exp (-K,°t) 

Expanding to meet the initial conditions 


ons) 
Up(Xy 0) ~ u(x) = =] = Bon cos (nex) 
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Using the orthogonality property of these eigenfunctions: 


fr (-1) cos (ns) dx = = a cos ( maps) ‘ 


cos (nmx) dx=> sin yo ery 
ii 
2 mit/ 2 > 


Eigenfunction constants: 
a. 2x «2 sin (m/2 

mrt/2 
Transient Solution: 


Un( x,t) = <2 sin (on/2) cos (nnx/2) exp [-(nn/2)*t | 


Total Solution to UA2: 





n(nn/2) cos (nix/2) exp [-(nio/2)*+ ] 


oo 
tha(x,t)iel -02 Sei 
: nel n/2 





L4.2 


L4.3.3 SOLUTIONS OF UB1 and UC1 
The total solutions to problems UB] and UCl are given in 

equations L4.3 and L4.4, respective:y. The development of the UB1 
solution is presented below; the UCi solution development is identical 
to UBl except for the obvious sign change. 
Statement of Problem UBl: Step L4.1 

Fe -wuee aa gC) 

(2,t) = U_, (t) 


u(o,t) = 0 


u = 0 


u(x,0) = Q 0 t 
ue= 0 


=52= 








Transient plus Steady State: Step L4.2 


u(x,t) = u(x) + Un(x,t) 
Osx YO aU = Uy = Une 


u.(o) = 0 uplo,t) = 0 Up(x,0) = -u,(x) 
u.(1) = 0 un(l,t) = 0 — up(x,00) = 0 
Steady State Solution: Step L4.3 


Transient 


u(x) = A+B exp (x/q?) 
u,(o) = 0 gives A+B=0 
u.(1) = 1 gives A+B exp (1/q?) =1 


solving the two A, B equations gives 


u(x) =l-e x/92 
1 = exp (1/q? 
Solution: Step L4.4 
Un(x,t) = X(x) + T(t) 
T(t) = exp (-K*t) 
X_ (x) = (A, cos Dx + A, sin Dx) exp (x/2q?) 
Uplo,t) = 0 gives A, = 0 
un(1,t) = 0 defines eigenvalues sin D=0 gives D. = nq 
‘ e a ji : - 2 
where: D2 = eee 1 = (nq) 


solving for the eigenvalues 
K 2@=1 + (ntq)* 
“ ies 

EKigenfunctions are: 
a, sin (nrx) exp (x/2a?) exp (-K,?t) 
Expanding te meet the initial condition: 
un(x,o) = -u,(x) = - [1 - exp (x/a?)7/c 

= $ a sin (nx) exp (x/2q?) 

n=] 7 


where: C =1 = exp (1/¢?) 
ahee = 





Using the orthogonality with weighting function property 


of the eigenfunctions: 


fr [exp (x/a?) ~ 1] exp (~x/ 2a) sin(oorx) dx 


oo 
Rey 2 ia sin (mmx) sin (nx) dx op 
Rae" O 


m 
-1) (mr) exp (#1/29?) = se 
[(mir)* + (1/2a 2 


Kigenfunction constants: 


a, = 2(=1 “enam) exp (-1/202 
nty)~ + (1/20 
Transient Solution: 


foe) Nn 
up(x,t) = 2 E, dol) (nr) exp [(x-1 202] sin (nrx) ° 
Li + (1/2a 


exp [-(1/4a? + a? n?n?) t] 


Total Solution to UBl: 


| u(x,t) = [l= ex an + 2 exp [-t/4a? + (x-1)/2a?]} 
1 = exp (1/a*) | 


mm ~] ast obs (anrr)*t }] sin (nx) 


L4, 3 


Total Solution to UCl: 


u(x,t) = [1 = exp (~x/g?) | + 2 exp [- t/t? + (1-x)/2a?]° 
i - exp (-1/a* 


co nN 
“~1) (nr) ex ~ (anrr)*t]} sin (nr) 





L4.4 
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L4.4 NUMERICAL CALCULATION AND GRAPHICAL PLOTS OF THE SOLUTIONS 

Bach of the solutions to the four problems of Table L4.2 contains 
an infinite series in eigenvalues and eigenfunctions. This section 
describes briefly a digital computer program which was written and 
utilized to numerically calculate the values of u(x,t) for q = 1.18 
at ten spatial points and nine time points. These values were 
determined for the purpose of making graphical plots of the solutions. 

The computer printouts and description of the program are given 
in Appendix A4, The computer output plots have been used to make the 
following graphical presentations of the four solutions which are 
presented in Figures L4,3 through L4.6. It is helpful to visualize 
these solution curves in terms of Figure L4.2. Each solution curve 
from left to right represents a spatial distribution at a given time, 
The curves can be seen to approach their steady state values as, for 
example, the linear distribution in UAl or the curved exponential 
distribution in UCl. 

The digital computer program in Appendix A4 computes all four 
of the solutions to three decimal place accuracy. It was found that 
at all spatial points and time points the number of eigenvalues and 
eigenfunctions, i.e. terms in the series, necessary to achieve three 
decimal place accuracy was in all cases less than ten (10). 
L4.5 SIMILARITIES BETWEEN CHAPTER L1 SOLUTIONS AND CHAPTER L4 SOLUTIONS 

It is interesting at this point to investigate the similarity 
between the standard analyticai solution methods of this chapter and 
the general solution technique presented for the LPCM in Chapter Ll. 
This relationship can best be seen by relating each Step in this 


chapter to the corresponding equations in Chapter Ll. 
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Figure L4,3 = SOLUTION UAL 
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Figure L4.4 - SOLUTION UA2 
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Figure L4.5 = SOLUTION UB1 
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Figure L4,6 = SOLUTION UC1l 
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The problem statement Step L4.1 or, more generally, Figure L4.1 
corresponds to Figure Ll.1, the statement of the LPCM. The transient 
and steady state conditions are similar to equation L1.15 and the 
boundary conditions L1.8 through L1.13. In both cases, the steady 
state solution, Step L4.3 or L1.5, L1.8, and L1.9, represents the 
solution to a two point boundary value problem. This is easily 
found in the simplified equations but requires a computer program 
(such as LBVP) for the general case. 

The transient solution Step L4.4, is, of course, the major step 
in solving the simplified examples as it is, also, in solving the LPCM 
in general. In the simplified examples, the eigenfunctions X_ (x) 
are expressible as closed-form functions with very obvious properties 
which can be used to satisfy the boundary conditions and to immediately 
solve for the eigenvalues Ke In the general LPCM, however, once the 
equation L1.25 has spatially varying coefficients, the eigenfunctions, 
in almost all cases, cannot be expressed in closed form but must be 
represented by infinite series. The infinite series representations 
of the eigenfunctions can be found by a number of techniques, such as 
the power series Method of Frobenius (H-12, Ch. 4), but these methods 
have several major disadvantages in this application. 

The first major disadvantage to infinite-series eigenfunction 
representations is that their properties in terms of satisfying 
boundary conditions are usually not obvious. Secondly, these methods 
lend themselves best to specific example solutions and become very 
difficult to talk about in terms of a general analytical solution 


technique because of the “special cases" involved. Finally, the 
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eigenfunction expansion to meet the initial condition is analytically 
difficult, in general, when these infinite-series eigenfunctions are 
involved. For these reasons, two other analytical techniques for 
finding the eigenvalues and eigenfunctions have been investigated. 

These two techniques are the Prufer Substitution (B-14, Ch. 11) 
and Fredholm II - Integral Equation theory (Appendix Al). The Prufer 
Substitution technique is a transformation applied to the second- 
order differential equation L1.25 which results in two first order 
ordinary differential equations whose solutions define the eigenvalues 
and eigenfunctions. The resulting first order equations are very 
useful in showing the properties, such as existence, ordinality, 
separation, orthogonality, etc., of the eigenvalues and eigenfunctions, 
but the differential equations are analytically untractable in terms 
of specific solutions and did not seem to lend themselves to consistent 
approximations. For these reasons the Prufer Substitution has been 
discarded as a possible analytical solution technique and the integral 
equation technique has been developed in this thesis. 

The final part of the transient solution Step L4.4 is the eigen- 
function expansion to meet the initial condition L1.12. The key to 
success in this expansion for the simplified cases is the orthogonality 
of the eigenfunctions with respect to a weighting function over the 
same interval as the boundary conditions are applicable. This 
orthogonality condition is expressed by equation L4.5. Of course, the 
sine and cosine functions of the simplified LPCM examples satisfy L4.5. 
The real utility of the eigenfunction orthogonality condition in the 
Simplified examples lies in the fact that the constants A. can be each 


explicitly expressed in terms of n. 
-161- 





b 


I. w(x) X (x) X, Cx)dx = § Cs L4.5 


mn 


where: C_ = constants # 0 
ban 70 for m#n 
Son =1 for m=n 

One may well wonder why the orthogonality condition was not 
utilized in solving for the constants C in equation L1.57, when this 
is one of the most important properties of these functions. The reason 
is that the functions X (x) are the result of transformations L1.44 
and L1.45 applied to the solutions W,(z) obtained from Appendix A2 
or Appendix A3 and, as such, are not easily expressed in general 
form, let alone integrated, For specific cases of the LPCM, if the 
functions X, (x) are expressible in a form easily evaluated by inte- 
gration, then the set of equations L1.59, only with the orthogonality 
constants instead of the functions evaluated at points — will be 
individually solvable because the matrix X will be diagonal. In that 
case, a better solution technique would be to use the orthogonality 
condition rather than the discrete-point approximate solution technique 
employed in equation L1.58. 

The fact that the functions X (x) satisfy the orthogonality 
condition, in the general case, is a result of the characteristics 
of the original eigenvalue problem Ll.25. If the transformation 
relations used in subsections Ll.5 and L1.6 satisfy continuity con- 
ditions and result in the regular Sturm Liouville problem L1.40, then 
the kernel (Al.17) of the integral equation is guaranteed to be con- 
tinuous. In addition, it has been shown by Lovitt (L-8, pp. 181-182) 


that for the special case of R(x) = 1 in equation L1.40, the kernel 
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Al.17 is symmetric and satisfies equation L4.6. The kernel symmetry 
K(z,s) = K(s,z) 14,6 


guarantees the existence and reality of the eigenvalues and the 
orthogonality of the eigenfunctions, Of course, in the simplified 
LPCM examples of this chapter the function R(x) is always one (1). 
In summary, this chapter presents analytical and graphical 
solutions to four simplified cases of the LPCM. The techniques for 
finding these solutions is then related to the general, integral- 
equation solution technique for the LPCM presented in Chapter Ll. 
The next Chapter L5 presents a brief discussion of several alterna- 


tive solution techniques which might be applied to the LPCM. 
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CHAPTER L5 

OTHER SOLUTION TECHNIQUES APPLICABLE TO THE LPCM 

This thesis presents the LPCM as a dynamic model of the composi- 
tion behavior of a binary plate distillation column and suggests an 
integral equation solution technique for it. The speed and precision 
of the integral equation solution technique has not been demonstrated 
completely in this thesis, however, Therefore, this chapter describes 
briefly the basic principles behind several other numerical methods 
for the solution of parabolic partial differential equations which 
might be applied to the LPCM for the purpose of comparing numerical 
solution accuracy and speed with the integral equation method. 

Numerical methods for the solution of ordinary and partial 
differential equations are extremely useful in the study of distilla- 
tion columns since both the discrete-plate and the continuous-spatial 
equations are impossible to solve analytically except in the simplest 
cases. Even when analytic methods are applicable, it is often the 
case that the solution is expressed in the form of an infinite series 
rather than in closed form. When this is the case, it may sometimes 
be less time consuming to apply a numerical method directly than to 
evaluate a series to some desired degree of accuracy at each point 
of interest, This was not true in the simplified cases of Chapter L4 
because each of the series converged so rapidly, i.e. less than ten 


terms e 


There are a large number of numerical methods which can be applied 


to any particular distillation column model, such as the LPCM. The 


choice of a particular one may depend on the form of the LPCM for any 
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particular column. ‘This chapter presents the basic principles behind 
the explicit numerical methods used for solving parabolic partial 
differential equations, applies these methods to the LPCM, and presents 


the resulting recursion relations, 


L5-1 NUMERICAL SOLUTION TECHNIQUES FOR PARABOLIC PARTIAL DIFFERENTIAL 
EQUATIONS 


Parabolic partial differential equations are usually solved by 
the use of step-by-step finite difference methods (B2.4 = Numerical 
Solution Techniques). This method consists essentially of defining 
a regular (usually rectangular) mesh and replacing the differential 
equation by a difference equation defined on the nodes of the mesh. 
Most of the variations in types of methods result from different 
derivations of the difference equation from the differential equation. 

In the distillation column models described in Section 2 (M), 
the solution sought is a function u(x,t) where both u and x are in 
the 0 to 1 interval and t is greater than zero. A possible mesh for 
use in equations of this type is presented in Figure L5.1. In this 
mesh the solution domain is divided into intervals of width h in x and 
1 int. The two increments h and 1 need not be the same and in some 
cases may be changed over different portions of the domains. The 
ratio of these two increments is constrained by stability requirements 
for any particular equation and method of solution. 

Once the mesh is set up the next step is to consider u(x,t) only 


at the mesh junction points as in equation L5.1. The parabolic LPCM 


mux, &) ~ u(ih, iL) |Z ies L5e1 
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Figure L5.1 - TIME-SPACKE MESH FOR NUMERICAL SOLUTIONS 





can be solved by either explicit or implicit methods. Explicit methods 
use only previous and present values of u, such as uw 3 Us jeyr 
to compute the future value Us jel° Implicit methods use previous, 
present, and future values of u to approximate the derivatives at the 
present point (ih, jl) which are then used to refine the present value 
of ue Both methods are applicable to the LPCM, but only explicit 
methods will be discussed here. It is expected that explicit methods 
would be adequate for solving most cases of the LPCM, but if not, 
then implicit methods (B2.4 - Numerical Solution Techniques) could be 
employed, Explicit methods are generally easy to implement and rapid 
to execute but are relatively unstable compared to implicit methods, 
Most finite difference approximations are derived from a Taylor 


Series expansion of u about the value u. Only as many terms in 


iyJ° 
the series are retained as necessary for an accurate solution. A 


Taylor Series expansion in t with only the first order terms retained 
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is presented in equation L5.2. Similarly, Taylor Series expansions 
in x with up to the second order terms retained are presented in 


equations L5.3 and L5.4. 











Be a 4,5 * oa 1 15.2 
ee 

u eu, 2 - du h + 02y h3/2 15.3 

i-1,3 1g dx 9x? 
1g OA 

ee, 02 2 
Mia 5 % 4,5 + a h “ta h?/2 15.4 
O* 155 1. 








L5-2 APPLICATION OF EXPLICIT METHODS TO THE LPCM 
The basic idea is to use stepwise expressions for u and its 
derivatives in the LPCM equation, here presented in expanded form is 
L5.5- In that case the above expansions must be solved for the 
P, (x) ~ + P3(x) a + P,(x)u = e 15.5 
approximated derivates of u in terms of the previous and present 
values of ue. These solutions will depend upon how many terms of the 


expansion are used. For the above expansions the approximated values 


of u and its derivatives are given by equations L5.6 through L5.9. 





su, . fone 
u ig 5 
eu -~ Dis) ae VieL. j Goer 
x ; 
1,J A 
02) = Us+1, 3 ~ au, Ps Us-1,3 b5eo 





gx* 


eae he 





du =u, .-u, ; L5.9 


Substituting these approximate relations into the LPCM equation and 
then solving for Us4y 5 gives the recursion relation L5.10 for the 
9 


LPCM, where the polynomial coefficients are evaluated at x = ih. 


Us47,5 7 (1 - A, - B, + Cus ~ Ayu, 


sly J 
L5.10 
~ CoUs ) jel 
Wheres A, = 1 - Ph/P, 
= 2 
Be ane /P, 
athe 
Crt a /(P,1) 


P, (x) are defined in equations Ll. 30 
through Ll. 33 
Using L5.10 and the previously derived equations a numerical iteration 
scheme can be set up for solving the LPCM along the mesh. 
L5-3 BRIEF SUMMARY OF SECTION 3(L) 

Section 3(L) concentrates on the Linear Polynomial-Coefficient 
Model (LPCM) which is derived in Section 2 (M) to represent the 
transient behavior of a binary plate distillation column. Chapter Ll 
presents an integral equation procedure for analytical solution of the 
LPCM in general using Appendices Al, A2, and A3, Chapter L2 applies 
the solution technique to a reduced form of the LPCM having one spatially 


varying coefficient and one constant coefficient. A suggested scheme 
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for a computer program PLPBV to evaluate the solution is presented 
using Appendices A5 and A6, Chapter L3 evaluates the LPCM for several 
cases from the steady-state representation of an example column, 
Chapter L4 solves analytically and evaluates numerically using Appen- 
dix A4, four simplified examples of the LPCM: two cases of the best 
equation and two cases of the Taylor Diffusion Model. The analytical 
solution steps for the simplified examples are then related to the 
analytical steps of the integral equation technique of Chapter Ll. 
Chapter L5 discusses briefly several possible numerical methods for 
solving the LPCM and derives the detailed recursion relations to 


solve the LPCM by an explicit numerical method, 
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SECTION 4 


SUMMARY AND CONCLUSION (S) 


Sl SUMMARY AND CONCLUSIONS 


S2 AREAS FOR FURTHER STUDY 


" THE PURPOSE OF COMPUTING IS INSIGHT, NOT NUMBERS" — R, W. HAMMING 
Grteals) 


THE CENTRAL PURPOSE OF THIS THESIS IS THE DEVELOPMENT, PRESENTA- 
TION, SUGGESTED SOLUTION TECHNIQUE, AND EVALUATION OF THE LINEAR 
POLYNOMIAL-COEFFICIENT MODEL (LPCM) OF THE DYNAMIC BEHAVIOR OF A 
BINARY PLATE DISTILLATION COLUMN. THIS SECTION SUMMARIZES THE MAJOR 
ARGUMENTS PRESENTED IN THIS THESIS RELATING TO THE LPCM AND EMPHASIZES 
A LARGE NUMBER OF AREAS FOR FURTHER STUDY USING THE LPCM AND THE 


INTEGRAL EQUATION SOLUTION TECHNIQUE PRESENTED IN THIS THESIS. 


SO 





CHAPTER Sl. 

SUMMARY AND CONCLUSIONS 

This thesis is summarized in the format of the Ghapter and 
Appendix Relationship Diagram presented in Chapter C2, 
Sl.l THE MAIN READING PATH OF CHAPTER C2 

Starting from the basic principles of distillation this thesis 
develops a discrete-plate dynamic model for the composition behavior 
of a binary plate distillation column. The fact that discrete models 
characteristically have large solution times leads to the search for 
a faster model and to the investigation of a continuous-spatial 
dynamic model derived by treating the plate number in the discrete 
model as a continuous variable. An investigation of possible solution 
techniques for the continuous-spatial model leads to the conclusion 
that for any hope of an analytical solution, the nonlinear continuous- 
spatial model must be linearized. The representation of the spatial 
coefficients of the linearized continuous model by general n-th degree 
polynomials defines the Linear Polynomial-Coefficient Model (LPCM). 


The basic equation of the LPCM is presented here as equation Sl.l. 


Q?_(P(x)u] +9. [P,(x)u] =u $1.1 
ox? : - a Z Ot 


Where: P, (x): and P, (x) are polynomials in x. 


Two-point boundary conditions and step initial 


conditions are specified for u and Uy 





The steady-state solution u(x) of the LPCM is obtained by solving the 


two-point boundary value problem. The transient solution is derived 


Sh ie 








by converting the spatial ordinary differential eigenvalue problem 
resulting from separation of variables to a Liouville Normal-Form 
equation and further transforming to a homogeneous Fredholm II integral 
equation. The solutions of the integral equation then define the 
eigenvalues K, and the eigenfunctions X, (x). An eigenfunction expan- 
sion to meet the initial condition then defines the eigenfunction 
constants C,. The total solution to the LPCM is then given by equation 


51.2. 


= : v.32 
u(x,t) u(x) + BO4X 0) exp(-K, *t) sil 





The integral equation solution technique used to develop equa- 
tion 51.2 is described in detail by applying it to a first-degree 
polynomial model (PLPBV) and by suggesting the type of computer 
program which must be used to evaluate the solution. The individual 
steps in the solution technique are described analytically by applying 
the transformation relations to PLPBV and are described numerically 
by a flowchart of the necessary computation steps and by a very basic 
Fortran IV digital computer subroutine showing some of the details of 
the computation steps. Preliminary tests with the PLPBV program 
suggest that it may be possible to generate complete column solutions 
in less than one minute for 20 spatial points and 9 time points. 

The complete model of a binary plate distillation column using 
the LPCM consists of two equations of the form of 51.1 and their corres- 
ponding two-point boundary conditions. The steps necessary to derive 
the LPCM for general distillation columns are presented. The details 


of these steps are then applied, fora second degree model, toa 
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5-plate and an ll-plate column and the resulting models are presented 
aS examples. 

Four simplified versions of the LPCM in the form of the heat 
equation and the Taylor Diffusion Model are defined as simplified 
versions of the LPCM. These models are then solved analytically, a 
computer program used to sum the series is developed, and the solutions 
are presented graphically and numerically. The standard analytical 
solution methods for the simplified models are then compared to the 
general solution technique for the LPCM. 

There are a large number of numerical methods which can be 
applied to any particular distillation column model. Several possible 
numerical methods for solving the LPCM are discussed, and the detailed 
recursion relations to solve the LPCM by an explicit numerical method 
are derived. 

The thesis is summarized briefly, and a discussion of a large 
number of areas for further study is presented. 

Sl1.2 THE AREAS OF OVERALL PERTINENCE OF CHAPTER C2 

A bibliography containing 352 references of which 202 pertain to 
distillation column dynamics and control is presented. These references 
are then related to four major areas pertinent to the thesis. Some 

1. General Tneory of Distillation 
2. Distillation Column Dynamics 
3. Distillation Column Control 
4, Mathematics and Computation 
general notes, descriptions, and comments pertaining to most of the 


references are included. 
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Several opinions on the philosophical aspects of observing 
reality as related to modeling a distillation column and on the desire 
to achieve profit as related to the cost of separation of components 


are presented. 
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CHAPTER S2 


AREAS FOR FURTHER STUDY 

This chapter presents,in the form of suggestions, a list of 
possible topics for further study or areas for further research in 
the order in which they are encountered by following the main reading 
path of Chapter C2, It must be emphasized that the LPCM and the 
integral equation solution technique are suggested by this thesis 
to be valuable tools for modeling distillation columns, but the vali- 
dation of this suggestion can only come from further study. 
S2.1 SECTION M 
1. Using the general discrete equations developed from considering 
all four (or part) mass balances per plate, develop continuous-spatial 
models and then LPCM models for the general binary plate column. 
2. investigate other methods of converting discrete models to continu- 
ous models which may give continuous models of greater accuracy. 
3. Investigate analog computer solution techniques for the CSE and 
the nonlinear polynomial-class CSE. 
h, Find in the literature or develop transformations applicable to 
the CSE.. 
5. Apply the Fixed Point Theorem to the CSE, 

ae Show that the CSE is a contraction mapping, 
be Investigate techniques for finding the fixed point for the CSE, 

6. Investigate solutes techniques for polynomial-class partial dif- 
ferential equations to see if nonlinear approximated CSE*s are more 


easily solved and less accurate than the general CS. 
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7- Investigate the validity and accuracy of models linearized about 
operating points other than the initial steady-state. 
S2.2 SECTION L 
8. Investigate the validity of the LPCM as a distillation column 
model using programs developed and based upon the integral equation 
solution technique of this thesis. 
ae Complete PLPBV including all special cases, 
b. Develop main modeling program to start from the equilibrium 

curve and the column physical characteristics and to end up 

With complete LPCM and solutions using PLPBV. 
9. Develop the digital computer programs for the second-degree (and 
higher) versions of the LPCM and use them to compute column solutions. 
10. Develop numerical solution program for the LPCM and compare 
solution time with the analytical solution programs. 
ll. Extend the step solution capability of the LPCM to approximate 
the response to any arbitrary input function of time, 
12. Investigate the relative magnitudes of the eigenvdlues for several 
variations of the LPCM in the interest of seeing how many are really 
necessary to completely characterize column behavior. 
13. Apply separable kernel solution procedure to the integral equa- 
tion resulting from the LPCM and evaluate the accuracy. 
14, Investigate in detail the properties of the function M(z) for 
several variations of the LPCM, 
15. Apply the integral equation solution technique to the two-point 
steady-state boundary-value problem and combine with the transient 


portion, 
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16. Develop an LPCM modeling technique for the case when the relative 
volatility is not constant. 

17. Develop an LPCM modeling technique for the case when the Murphree 
plate efficiencies vary and are not necessarily unity. 

18. Investigate in detail the symmetry and continuity properties of 
the integral equation kernel K(z,s) for various functions M(z). 

19. Investigate implicit numerical methods for solving the LPCM and 
compare solution times to the analytical methods. 

20. Investigate the possibilities for partially analytical - partially 
numerical solutions to the LPCM (hybrid solutions). 

$2.3 EXTENSIONS 

21.- Develop distiliation column control schemes particularly suited 
to using the LPCM. 

22. Investigate the possibilities for partially discrete - partially 
continuous models for distillation columns (hybrid models). 

23. Extend the application of the LPCM techniques to packed columns. 
24, Extend the application of the LPCM to continuous models of multi- 
component distillation columns. 

25. Combine spatial LPCM and time LPCM solution techniques to solve 
partial differential two-point boundary-value problems in which the 
partial differential equation is second order in both time and space. 
The result of the application of separation of variables will be a 


spatial LPCM and a time LPCM. 
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SECTION 5 


APPENDICES (A) 


Al 


A2 


A3 


AY 


A5 
A6 


CONVERSION OF A LIOUVILLE NORMAL-FORM EQUATION TO A 
FREDHOLM II =< INTEGRAL EQUATION 


SOLUTION OF A HOMOGENEOUS FREDHOLM II - INTEGRAL EQUATION 
WITH A SEPARABLE KERNEL 


SOLUTION OF A HOMOGENEOUS FREDHOLM II = INTEGRAL EQUATION 
BY CONVERSION TO LINEAR ALGEBRAIC EQUATIONS 


DIGITAL COMPUTER PROGRAM AND OUTPUT FOR EVALUATION OF 
PROBLEMS IN CHAPTER L4 


SUBROUTINES USEFUL FOR COMPUTING LPCM STEADY-STATE SOLUTIONS 


SUBROUTINE PLPBV 


“MATHEMATICS POSSESSES NOT ONLY TRUTH, BUT SUPREME BEAUTY — 


A BEAUTY COLD AND AUSTER®, Likt' THAT OF SCULPTURE, * © © SUBLIMELY 


PURE, AND CAPABLE OF A STERN PERFECTION SUCH AS ONLY THE GREATEST 


ART CAN SHOW." — BERTRAND RUSSELL 
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APPENDIX Al 


CONVERSION OF A LIOUVILLE NORMAL-FORM EQUATION TO A 
FREDHOLM II ~ INTEGRAL EQUATION 


This appendix presents the intermediate steps necessary to the 
conversion of the Liouville Normal-Form equation (B-14) and its 
boundary conditions to a Fredholm II - Integral Equation (L-8,H-13). 
The conversion proceeds via integration of equation Al.1 and 
application of boundary conditions Al.2, Al.3, and Al.4 to integral 


equation of the form of Al.5. 


a - [M(z)-)]W = 0 Al.1 
D, 4 W(c) + DjoW(c) = 0 Al.2 
Do, W(o) + DopW(o) = 0 A1.3 
DizDe2 - Dy2Dz # 0 i 
W(z) = fo K(z,s)W(s)ds Al. 5 


Integrating Al.1 from 0 to Z twice and making a change of 


variable in one of the integrals yields Al.6. Application of the 


W(2) = f% [M(2-s)-a]W(s)as + W(o)z 


+ W(o) A1.6 
second boundary condition Al.3 to Al.6 then gives Al.7. The 


W(z) = i. [M(z-s)-) ]W(s)ds + [1-Dp,2]W(o) wi. 7 


Doo 


remainder of the conversion procedure consists of determining 


Ne 





the constant W(0) and placing the resulting equations in integral 
equation form. 
The first step in evaluating W(0) is to integrate equation 


Al.l from Z to C, giving equation Al.8. Evaluation of Al.8 at Z=z0 


Wie)xcew(n) + S° [M(2)-nJw(2)ae A1.8 


and of Al.7 at Z =C gives Al.9 and Al.10. Placing the boundary 
fe) ° Cc 
W(c) = Wo) + f [M(z) =, |W(z)dz Al.9 


Cc 
W(c) = JjLu(e-s) -a Jw(s)ds +{1 - Dor ¢ Jw(o) A1.10 
22 
conditions in the form of Al.11 and Al.12, substituting them into 
Al.9 and Al.10, and equating them gives Al.13, an equation solely 


in terms of W(0). 


D 0 

W(c) = = 12 W(c) Al.11 
Diy 

0) D 

W(o) = = _21 W(o) A1.12 
Doo 


J [a (e-s)-n]w(s)ds +{[1 - Por © jw(o) 
22 


~ - Pra [- Por wo) + JO[m(z)-aJw(z)az] 1.23 
be 


Solving for W(0) gives Al.14 and Al.15. 


W(o) = 2 J°[(M(c-s) + 2 m(s)) - 422), 
e “ih "1 


W(s)ds Al.14 
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Gham clad, + “gal 14.0 Al.15 
Di Foe, Se 
If equation Al.14 for W(0) is now substituted back into equation 
Al.7, then equation Al.16 is an integral equation in W(z). This 


equation Al.16, however, is not yet in the desired form of Al.5. 


W(z) = f7 [M(z-s)-)] W(s)ds + 1_ (1-21 z) 
o cS) as 


so [(w(c-s) + 722 M(s)) - (e212),] 
sigh Hh 
W(s)ds Al.16 


If the kernel function K(z,s) is expressed as the sum in equations 
Al.17, Al.18, Al.19, and Al.20,then the conversion to the integral 


equation is complete. 


K(z,s) = K, (2,8) + K,(2,s) 
Sku S>Z Al.17 


K, (z,s) = M(z-s) + 1 (Doq*Do z) 
1 


(D, ,M(c-s)+C, oM(s)) 


-, [1 + 1. (D,5-Do12) (Dy +D, 9) J Al.18 
1 


K5(z,s) = a (9 ,07Pa) [D, M(e-s )+D, .M(s) 
if 
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Thus, the complete integral equation representation of the 
Liouville Normal-Form equation is summarized in equations A1l,20 


through Al.24. 


W(z) = ie K(z,s)W(s)ds Al.21 
K(z,s) = K, (zs) + K,(z,s) Al. 22 
K, (z,s) = M(z-s) - ) + K,(z,s) s<Z Al.23 


K5(z,s) = g Pz2rPaa*) LD, ,M(c-8)+D, oM(s) 


= (D,,+D,5)aJ Ss > Z Al.24 


Now, the kernel given in Al.22 must have two properties in order 
for solution procedures to be applied. It must be continuous at s = z 
and it must be symmetric. The continuity aspect implies that A1l.23 


can now be written as Al.25 and Al.26. In this case a new function 
K, (zs) = M(z-s) - M(O) + K,(z,s) Al.25 
K, (2,8) = \K4(z,s) Al. 26 


K(z,s) has been introduced, where K,(z,s) is given by Al.27 and 


Al.28. Utilizing once again the continuity requirement and the 
K5(z,s) = -M(z,s) + M(O) + AK(z, 8) Al.27 
K, (2,8) = Ky(z,s) A1,28 


definition of Ky, (2,8) given in Al.28, then K), (Zs) can be written 
in terms of K,(2,s) as in equation Al.29. The remaining kernel 
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K, (2,8) = 1 = M(zes) + K3(z,s) Al.29 
M(O 


Where: M(0) # 0 


function K(z,s) can then be expressed as in equation Al. 30, and 


a general test for symmetry can be given by equation Al. 3l. 


K(z,s) = M(z,s) -1+1 (Doo - Do42) ° 
M(O a 


[D, ,M(e=s) + DjoM(s) - (D,, + D,>)}] Al.30 
Sap 


C) [M(z-s) - M(O)] + (Doo - D512) [Dy M(c-s) + Dy oM(s) 
- M(0) (Dy4 + Dy) J es (Do, - Ds) [D,  M(c=z) 
+ D, oM(z) - (D,, + D, >) M(O) | ac al 


The summarized equations for the integral equation representation 


are presented in Figure Al.1l. 


W(z) = oe K(z,s)W(s)ds 
K(z,s) = K,(z,s) + K) (z,s) 


S< Zz S > Zz 


K,(z,s) ™1 = M(z,s) + K.(2,s) S>2Z 
k uO 3 


K.(2,s) = eal 1+ 1_ (Do, = 05,2) 
, MCO or a 


-Prit(e-s) + Diot(s) - (yy + D,p)4 8 <2 
7"  “)y = = 


Figure Al.l = THE FREDHOLM II = HOMOGENEOUS INTEGRAL 
EQUATION 





-18 3- 








APPENDIX A2 
SOLUTION OF A HOMOGENEOUS FREDHOLM II- INTEGRAL EQUATION WITH 
A SEPARABLE KERNEL 
This appendix presents the developments necessary for finding 

the solution to a separable kernel integral equation, The integral 
equation A2.1 and separable kernel A2.2 are combined and a matrix 
solution technique is developed. The technique to be presented is 
standard throughout the literature when it is recognized that the 
following terms are equivalent to “separable kernel" A2.2, 

1. Separable Kernel (H=-13) (L-8) 

2. Degenerate Kernel (P-4) (G7) 

3. Kernel of Finite Rank (S-11) 

4, Pincherle-Goursat Kernels (T-3) 


5, Riesz-Schauder Equations (M-6) 

W(z) =} ie K(z,s)W(s)ds A2.1 
m 

K(248) =,2,44(2)®, (s) N22 


If A2.2 is substituted into 42.1 and the integral and summation 
signs interchanged, then equation A2.3 results, It can be seen that 
the integrated terms depend only upon S and are therefore constants 


A2.4 after integration. Thus, the solution to A2,1 can be expressed 
Tm c 
W(z) = A, 2s 24 (2) ni b, (s)W(s)ds A253 
C 
t~ 4 A2.4 
C, pe b, (s)W(s)ds 


in terms of the constants C, as equation A2.5, and the problem now 


becomes that of finding the Cie 
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W(z) = dB ay (2), A2.5 


The first step in finding the C, is the substitution of A2.5 


aL 
back into A2.4 giving A2.6. The terms involving the integral of 


4 ae equation A2.7, 


which can be calculated from the given kernel function A2.2,. The 


a, (s) and b, (s) are then recognized as constants A 


resulting equation A2.8 is a system of m simultaneous equations in 


m S 
Cy = A 2A 5 ae b, (s)a,(s)ds A2.6 
Ay; = I. b, (s)a (sds A2.7 


the m unknown constants C,. 
Equation A2.8 now represents a matrix eigenvalue problem, This 
Can be seen by expanding A2.8 and writing it in the form of A2.9, 


where 06 = 1/\. Equation A2.9 can then be written in the form of 


© 


oil 


Oe e8 8 CE 





GO#l A2Z.10 


A2.11, which is a matrix eigenvalue problem. 


(oI - A)C = 0 A2.11 
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The matrix eigenvalue problem A2.11 must now be solved. The 
first requirement for solutions to exist is that the determinant of 
the matrix be zero as in A2.12. This specifies an nee order 
polynomial in oO which when solved gives the m = real roots or 
eigenvalues G53 i =il1,m. For each of these Oy there must exist an 


eigenvector Ey such that A2.13 is satisfied. 


det (cI-A) = 0 A2.12 


(0, I-A)E, = 0 A2.13 


Once these eigenvalues and eigenvectors are known, the first 
solution can be expressed by A2.14., It must be emphasized that there 


are m - solutions W, (2) to this problem, one for each eigenvalue. 


m Io 
af 


Thus, the eigenfunction solutions to the equation Az2.1 are given 


an A2Z.15. 


m 
W, (z) = 2 4585(2)8, A2.15 
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APPENDIX A3 


SOLUTION OF A HOMOGENEOUS FREDHOLM II-INTEGRAL EQUATION BY CONVERSION 
TO LINEAR ALGEBRAIC EQUATIONS 


This appendix presents an approximate solution technique for 
determining the eigenvalues and eigenfunctions of an integral equation 
A3.1 by dividing the integration interval into equally spaced in- 


crements and changing the integral to an algebraic summation. The 


W(z) = J K(2,8)W(s)as A3.1 
integral equation eigenvalue problem is thus changed to a matrix 
eigenvalue problem. The approximate solution W(z) can then be 
expressed in terms of the eigenvalues and eigenvectors of the matrix 
eigenvalue problem. References for this material are: (H-13) (L-8) 
(P-4) (G-7) (S-11) (T-3) (V-4) (S-14) (D-5) and (M-6). 
The first step in this solution technique is to divide the 


interval [o,c | into n equal parts of length § as in equation A3.2. 


6 = Boos 


a hi 


Next, the functions W(z) and K(z,s) can be evaluated along the 


interval as A3.5 and A3.6 at each of the points given by A3.3 and A3.4. 


K(16,56) = Ky, (i,5 = 2, 2, ves n) A3.3 
W(46) = W, (4 = 1, 2, ¢**, n) A3.4 
Z, = 16 (4 = 1, 2, °**, n) NERS 
Ss; = 38 (j =1, 2, °°, n) A3.6 


With these relationships established, the integral in A}3.1 is 
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approximated by a summation over the interval [0,c] represented by 


the indices [1,n] as in equation A3.7. 
W : K, .W A 
ey ie ran om 


Equation A3.7 represents a matrix eigenvalue problem similar 
to equation A2.8 in Appendix A2. This equation A3.7 can be placed 
in eigenvalue form A3.9 by defining 6 as in A3.8, similar to equations 


A2.9, A2.10, and A2.11 in Appendix A2, 


ee A3.8 
r§ 
(Q@I-A)W = 0 A3.9 


Once the eigenvalues O, and eigenvectors Ey for equation A3.9 
are found, the approximate solution to equation A3.1 can be expressed 


as A3. TOS 


W3(z) = Ago 2 K(2, 58) 3, , A3.10 
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APPENDIX A4 
DIGITAL COMPUTER PROGRAM AND OUTPUT FOR EVALUATION OF PROBLEMS IN 
CHAPTER L4 

This appendix presents a listing of the Fortran IV statements 
in the program and subroutines used to calculate the numerical values 
of the four analytical solutions presented in Chapter L4, The numeri-~ 
cal results are first listed, then plotted, and finally printed out in 
matrix form. Tests made with this program showed that in all cases 
less than ten terms in the series were used for three-decimal-place 
accuracy. 

There are three subroutines used in the program: MXOUT, PLOT, 
and LOC. Subroutine LOC was taken, less comment cards, directly from 
reference (I-3), the Scientific Subroutine Package. Subroutines MXOUT 
and PLOT, used here, represent significantly modified versions of the 
MXOUT and PLOT subroutines described in reference (I-3). The main 
program and subroutines are written in Fortran IV, and all runs were 
compiled and executed using the WATFOR compiler on the IBM 360/65 
computer at the M.I.T. Computation Center, The entire program 


compiles in 5.8 seconds and executes in 5.65 seconds. 
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C PRCGRAM TO CALCULATE 4 TRANSTENT RESPONSE S,UAL,UA2~,UBl UC 
C PROGRAM TO CALCULATE 4 ANALYTIC SOLUTIONS AS A TEST 
DIMENSION UA( 100) ,UAT( 100) ,UB( 100),UC(100) 


— mew cw cw eee eee eee ee ee eee eee eee ee — er ae oe o_o eee ie 


I eee ee ee eee ee eee ee = 


1 ay T=",F7.3) 


103 FORMAT(1X,*UAOKNE NAOT ONE PERCENT ACCURATE ,X=' ,F7. 3, 
| ad =*,F7/.3) 


A ce ee A a tl eee ne ee ee 


-— cr rr crm rr rw mc cr rm wm mm cmc mmm mcm cc ee es 


PI=3.14159265 
DELTT=0.025 
DELTX=0.1 
Memeeearrises 
BN=0.25*8 
WRITE(6,100) eee 
C TEIS SECTION ITS TO SET UP QUTPUT PLOTS AT OQ.1 INTERVALS 


i tt te el eel 


mem cre wr crc rr rr ww rc ww wr re wr we ei 


XP=XP+DELTX 
UAT) =XP 
UAT (I )=XP 


ems mm ww wre rr mm me me wes mm we we ww ese 


cere mec cs ee ee ee we we ee me we see es 





DO 10 J=l 29 


ee mmm rr rw ee ee ees esse a aa eee 


emer cm cr rm re ws we we ee i Ss Ses a es  - ee — ss se 


X=X*#DELTX 
NOw SET 


eee eee OK TD TS *NOweS Et 


C STEADY STATE CALCULATIONS 


tee ear cee ee ee ee ee eee ee ee ee eee ee eee eee nw ee se ss 


—— eee ee eee eee ee eee Oe eS SS Se 


UBOST=(1.0-EXP(AC#X))/1 .O-EXPLAC)) 

UCCST=(1.0-EXP{-—AC*¥X))/(1.0-EXP{—-AC)} 
C TRANSIENT CALCULATICNS 

SN=-1 Fae) 

UACTR=0.0 
Me _UATIR=0.0 _—_ 

UBCTR=0.0 

ATB=2.0*EXP({-0 -25%T+025*X-0.5) *AC) 


— ei ae ee —_——_— = — . SS ee ee 


eee eee 
me ei aS ess 


“1905 > se eee | a. 


mw we a es Ss ee ee 





rm rr ct eee —— —— ama oe wee ee oe — — a we 
— -_ 


eee eee ee —= 
— —_ ——————— mm mm ee ee we we ee ee ee 


DO 30 N=1,4C0 
IF (CIF IN-3)53,54,54 
53 RN=FLCAT(N) 


—_——s cee sc cm cw cw we ee es ee es ws ee es es we es es wees we ee ee ee eo ee i ee ee ee ae 


—_—_ own ie cm cr cr cr rm me ww we we ws es ee ee 


41 RTN=FLOAT (2*N-1 ) 
ATN=RIN*PI ~ 
COSN=COSCATN*X¥*0.5) 
EXPBN=EXP (~-RIN¥RTN*¥T XBN) 


—_—_—s—— sen cr cr cr crm rc mc wm mm ee ee ee 


— oe mm crm rm cr es ws es es  - 


IF (ABS(UATTR)—MIN) 44544,61 
61 TFC(ABS((UATTR-UATT)/UATTR)—-0.01)43743,42 
44 UATTR=0.0 
TEINETFING, 0 0,0, 2... OOO 
Meme CO TINAESINGANEX) OU 
EX PB=EXP{-RN*¥RN*XT¥#B) 
1F(1AG)45545946 
45 WAO=UAOTR 
______UAOTR=UAGTR+(SN*SINN*EXPBI/AN 
IF{ABSCUAGCTR) —MIN)48,48,62 
____ 62 TFUABS((UAOTR-UAC)/UAQTRI-0.601947947,46 
48 UAGTR=0.0 
47 1AQ=1 OTs rene 
IFIN=IFINGL 
46 IF(1B0)49,49,5C 


ee mm mcm rc mr mm ccm re ei ee 


ee anne gO ee I cl i 


UBOTR=UBOTR+{SN*AN*SINN¥EXPB) /SSUMSQ 
TFCABS(UBOTR )—-MIN) 521,52, 63 
63 IF CABS{((UBOTR-URO)/UBOTR I-0.01)51,51,50 


= = me ere eee ee 


Ss > - scm cm crm mr eee 





re eet Ss ee ND ee ee ae eee 


56 ITF(TA0Q)57,57,58 
fete Os lO3)X,t ee 


ee eee a) ee Pw ND LN Le 


So e1F(1B0)59.559,54 
59 WRITE(6,104)X%,_ 7% 
54 ATR=2 .0*UAQTR 
MONE SUAOSTHATR on ee 


o— rr ted 
creme ei ee ee ee 


ATTR=2 .O*UATTR 
Wee OeUATSTTATIR  — _.LKUL WLW Eee ee 


meme ee ee ee i ees ea Ser ee _ OO Ee I  — ee OO 


BTR=ATB*UBCTR 


MEMEEGTRSATC#HUBOTR 
UBONE=UBOST+BTR 


Meee UCOSTZECTR cc eee 


— <i 
i iil 


C SAVE VALUES FOR PLCTTING 
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——. ee cee ce er cm ce ee ee ee ee ee ee eee eee 
— om ce ee ec ce crc cr ccm cw eee i 


UA(I) =UAONE 
ee SS eTS=UATWOO LU 
WT =UBONG 000 IE EE SS TS 
UC (I) =UCONE 
I=1{+1 
SemeennrINue © 0 S55 
Zz 10 CONTINUE 
MeM@IP COMMANDS ..- © 0 SS 
C MATRIX PRINTOUT 


NQ=1 


—_— oe cc cr eo re ee cm mm cr mcm mm cm mr ww cm ew ee a eee 


NN eee 


mm cer rm es ce cc cc mp cc cc cm wer ec eee ei ee 


CALL PLOT{NO, UA,NP,MP,NL,NS 


smc cts ce ce ee eee eee ee ie oe 


mcm mr cr cr me cm cm eee eee i oe 


CALL MXOUT(NO, UAsNPyMP,MS,LINS, IPOS ISP) 
NC=NO+1 

CALL PLCTINO,UAT»NPyMP,NL,NS) 

CALL MXOUT(NO,UAT »NP,MP,MSeLINS, IPOS, ISP 


mr mm rm mm rr cr ce ee ee eee ee eee eee eee eee 


ecm ccc cre rr er crm rm cr cr we ee ee eee eee eee 


CALL MXCUT(NO, UB sNPyMPyMSeLINS, IPOSe ISP) 
NO=NO+1 

CALL PLCTINO, UC,NP,MP,NL,NS) 

CALL MXOUT(NG, UC,sNP,MPyMSoL INS, IPOS, ISP) 


esc ee mm we ce ww ww ee 


em rm mms mm crm emer cae eee a ees es 





ees eee omer cee ei le i ce ee ee ti me ee ee ee ee ee ee ee ee ee ee ee ee ee eer — eT ss = 2. 2 — CC = -_  _.CT 


ee eee eee ee ce es es es es cme re ce i re es ee i ee ie icc a ss en se eT ee Or ee 


ee eee es ee ee es ee ee a a ee es ee ee ee ee ee ae Se eS eS — eskK— SS ES = wa ee ee sh ,,T™TLTTD 


mmm eases eee eee 
emma se es 





SLBROUTINE PLOTINO »AyNyMyNLyNS) g 
c Me ok BR a MICHAEL Ne HAYES tak XR oe ok KK 


— 
_—_— ——_— ee ESE Eien rec ee me mm ee ee eee ee ee ee se ee eee 


DIMENSION OUT(101),YPR(11),_,IANG(9},ACL) 
PATEGER TDUM/*1*/,TANG/* 1% ,%2* 59 3%5'4*, 5%, 96%, *7', 


1*8*,#9#/ 
Oe i) 
C THESE LIMITED FORMATS ARE FOR 60 SPACE PRINTCUTS FOR 
€ THESIS USE. 


SE ——EEe————ESe See eee ee cm cm crm rc ce cr we ee ewe em we me ee eee wee wee ee 


1 FORMAT(1H1,28X,7H CHART 4135//) 
2 FORMATCLIH ,F8.391X o* ** g51LAL, ***) 
3 FORMAT(1H ) 

? FORMAT(1H ,10X ,36h* * * x * * as % 9 


rr —_— Se rc cm mc crm ec ccc cre me ee ee ee ee ee eee 


_— oe ee ee cm rrr mm mm mr wee ee eee 


C NTH ITS NUMBER CF SPACES TO BE USED, EITHER 101 OR 51 
NTH=51 
NLL=NL 
IF (NS) 16,16,10 


ce mets reece mmc mmc mm mew crm mre ce ccm ee ee ie = 


mmr cece cc cm ec wr cme me me em w+ 


TFC ACIS )-ACT0914514,11 
11 L=I-N 

LiL=J-N 

DO 12 K=1,M 


tet mmm mm mm em rt es ce ce ce ret ee + 


meme cmc em we mmm mem mw mr mmc wr cc tm ce cc ce eee eee 


A(LI=A(LL) 
12 A(LLI=F 
14 CONTINUE 


em mm rr se mmm rm ee ee asses eee 


mee ea se el aww ecm me ee ae 


18 NLL=50 
20 WRITE(65135NO0 
WRITE(6,7) 


mem cmc mr cm eee EE Er ee —_— 





mmm wc rr re wr te wc rm gs we wi es SO Se SS CO Oe? 


YMAX=—1.0E75 
YMIN=1.0€75 
M2=MEN 

— MeO S-41,M2. . eee 
IF (A(JS} «GT. YMAX) YMAX=A(J) 
IF (ACJ) obT.s YMIN) YMIN=AQJ) 


em ee ee ee eee eee ss 


40 CONTINUE 


YSCAL={YMAX-YMIN)/50.0 
XB=Ai1) 





rm er ee ee ee ee eee ee SS SS ES ee S= eI 


XPR=XB+F¥*XSCAL 
IF (XPR~-A{L1 150 750%70 


—— eee eee 
rete ee ee ee ee ee ee ee eee eee 





a TS eo eS cane oe ee ee ee ee 


———— een = _— en cc we ce eee ee eee ee eee ee 


60 


——_— ee rm mi cre rr eee ee ee ee ee 


I 


70 
8C 


eee ee ee ee ee 


ea ee 


ee ie 


ae 


me ee ee ee 


ee eel 


— ee ee eee ee 


comme ee -- 


JP=((AC(LLI-YMIN) /YSCAL)#1.0 
OUTCJPI=IANG(S) 
CONTINUE 


—_——_—_— rc mm cc mm ee ce eee ee eee 


————_—r—XK- Oe rr mm cr ww ss ee ee ee ee ee ee 


WRITE(6 3) 
T=f+1 


—— se cr cr crm cr crm cme crm ec ec rm ee eee ee ee ee ee ee 


—————————— crc ce rc cc rm rc crm ce cc cs ce ee ee ee ee ee 


WRITE (6,7) 
YPR(1)=YMIN 
DO 90 KN=1,9 


—_—_—_———— eee ec cr rr cr cr crc cr cm cr me cm rr ee ee ee ee Oe 


ee es ee ——— oe 


RETURN 
ENC 


eS cm me ccm errr mmm mcr cr cc cr cr cr cr ce ee eee eee ce 


em cr crm rm crm rr cr cr cr cr cs css cr cr cr cry crm ce eee ee 


rE ——_ 


ee de ee ee ee 


ce ccc ce cr eee 


meme ew eee 


mee, cmc iS ce rc mee cc ee ee 


te el 
mmm es eee 





! SUBROUTINE MXQUT (TCODEsAsNsMeMSyLINSy IPOS, ISP) 


o_o em mmr mc ee ew we we we ew ee wee ee ee ee 


DIMENSION AC1),B(8) 
1) FORMAT(1HO,5X, 7HMATRIX ,15,6X,13,5H ROWS, 6X,13, 
18H CQLUMNS,/) 
__ Zz é FORMAT(12X,8HCOLUMN 47(3X,13,10X)) = 
3 FORMAT(1H ) 
me! Hommel » (XeS4HROW »l3s7(E16.6)) == = = _ 
5 FORMAT(LHO, 7X,4HROW , 13, 7(E15.6)) 
6 FCRMAT{I1H1) 
T FORMAT(16X,13HSTORAGE MODE ,Il+s?X,7HGROUP ,12,/) 


WRITE (6,6) 


I I I a re th te a A eee 


mm tcc cmc mm em mm cr cm cc mc ee ee eee ee ee ieee 


LEND=(LINS/ISP)-2 
IPAGE=1 
10 LSTRT=1 
Secor EVOs1)ECODE,N,M ee 
WRITE(657)MSyI PAGE 
JNT=J+NEND-1 


i ele tt teil 


IPAGE=IPAGE+1 
31 ITE INT-M) 33 ,33 532 
32 JNT=M 
MS PCONTINGE = ee 
WRITE(68,2)(5JCUR,JCUR=J,JNT) 
______ IF ATSP-1) 35535940 
35 WRITE(6,3) 
40 LTEND=LSTRTI*+L END-1 
DO 80 L=LSTRTeL TEND 


DC 55 K=1,NEND 


i i i nite 


ema ee eee ee ee eee ee ee ee 


meee ee ee eee eee ee eee eee eee ee 


— ee ee ee ee ee ee er ees 


CALL LOC(LsJTsITINT sNoMyMS) 
B(K)=0 .0 


[FCI INT)50,50,45 
PARSING) UU a eee 


sm mmm cm me ee 


eed cm ee ee ee wea ae 


55 CONTINUE 


60 TF(TSP-1)65,65,570 


65 WRITE 6e4)L (BC UW), JWH=1,) KK) 


some cm mee eee ee we ee ses a ss ee ee ss 


tie 


80 CONTINUE 


GO 0 20 


emg 90~995,99 2 eee eee 


“95 RETURN —— 


END 











_- SUBROUTINE LOC(CLyJSsTReNeMsMS) 


c Re eB ie ee ae eae a ik eK Kae oe aie te aa ak ke Se a a ie ae ak ake ake ate at ake fe ake ak oe ae ae te age ake te ok ae 
es — 
JX=J 
eee EF (MS~1) 10.20.30 ~.. LL ae 
10 IRX=N*(JX-1)41X 
. SOOO. ee 
20 TFCIX-JSX) 22924524 
= + aK = 
GO TQ 36 
ee moe R X= Ke IGETX—1X)/ 255... Le eee 
GO TO 36 
= PRAT O ee ee ee 
be ( 1 5X) 36532,36 
= [x 
36 ITR=IRX 
ree WRN 8 ee, eee 
END 


em cm mm mmm meee ee ee SE s-2«s = 


meee ee 


mmm ms mmm cm mr rs we mm em mm eee ne 


ee rm mm ms ms me ty we cm wm ce mm crm crm, wy mm we mw ew ws ce ww ee mi me em as 


em mm ess a ema ee ee 


i eee eed 


ce mmm cme we mm a ai 


TL LL 


ec em wm wr wes me mw ww i ee EE ESsSeG S= ss ss eee er OO Oe ee 


ee ee 


— ee ee —_—— mm ee ee ee 
ee ie titi ined — —— — —_ — oom 


ee et es ee ee ee ee ee ee ee ee ee ee ee 00 SS 2 22:20 2 SS SS | Cee ee ee ee ee 


om i i lettre 
ee mm we ee ee es ee ee eee ee i eee a i ss ses ee eee ees 


—_ = stem ees wih meme em eee ee ee ee ee 
— eer ee ee ee ee mee ey ee ee ee ess es ee Lf ° i =} 


mm ee 
cows em em ee ee ee ee ee ee ee ee ee ee ee ee ee i se eee es oe =e er 


me ee 


em ee eee eee 





cm cr te crt cc cm eee ee eee ee ee eee eee 


<tc ce cr cr rt crm ec ee eee ee eee ee 


cr rr ret ct cc crm ee eee ee eee eee ee 


0.30595E-03 
0. 13734E-02 
0.50529E-02 
QO. 16391E-01 
0.32329E-01 
0.-10959E 00 
0.22939E O00 
0.4237CE 00 


0.37324E-03 
0.3569 1E-02 
0.4980 1E-02 
Q0.166C5E-0O1 
O-10548E 00 
O.22987E OC 
C.42378E 00 


O.21213-03 
0.993C1E-03 


Om BE37E-C2 


—— ce cc ec ee eee eee ee 


9-38026E-G2 
Q0.12842E-01 


O ub 8S E—-OZ 


emcee eee 


0.-26402E-01 
0.93069E-0O1 
0.20286E 00 
Q0.39023E O00 


0.39390E-01 
O.12817E 00 
O0.25789E 00 


om ccc ct cr ct ct ct cr cc rr rrr wt mec mm ee eee ee 


0.68921E 00 
O-1COQOE O01 


0.6051 GE 206 
O-10000E O01 


0.66130E 00 
CeLOOQOOE Q1 


O0.71638E CO 
Q-10000E Ol 


tet em etre cr crt tr ctr ce mm mm rer ccc cm cmc ee ee eee ie 


OF 9037 FE— O02 
0.22954E-01 
O.4/7477E-0O1 
0.89609E-01 


CslZT26E=01 
O23 2439 05-01 
0.48038E-O1 
O.8976LE-01 


0-.62542E-02 
0.16539E-O1 
0.35616E-O1 
0.69999E-01 


0.12850E-01 
Q.3136/7E-0O1 
0.62355e-O01 
O-11313E QO 


cme cr rr cm cr cr cc re cr crm cry ct tc cc cc cr mc te eee 


0.15724E O00 
0.2578C0E O00 


Os IZIZTE OG 
0.21850E O00 


Ce18987E CO 
O.30091E CO 


me cm cc ccc ree crc ce mmc cc tm mc mc ee eee 


O.57160E OO 
O.77730E O00 


0.39612E 00 
QO.5716CE 00 
O-77729E O00 
O.-1000CE O1 


0.34967E O00 
0.-52579E 00 
0.74534E 00 
QO-LOOQOOOE Ol 


0.44453E CO 
O-61703E CO 
C.80742E O00 
QO.1O000E O01 


em cr tm mm we ee ae ie i 


0.26589F-01 
0.59082E-01l 


0.48732E-O1 
0./6247E-Ol 


0.18346E-C01 
0.42445E-01 


0.37694E-O1 
Q0-80500E-01 


ems mcm mm ee mt mcm cc cc ee ee ee eee 


0.-10328E 00 
0. 16462E 00 
0.24762E 00 
Q0.35538E 00 


O.1080SE Q0 
O-167C7E OO 
0.24873E O00 
0O.35581E 00 


O./77256E-01 
QO.12827E 00 
O.20100E CO 
0.30057E 00 


QO.20731E 00 
O0.29986E OC 


me cmp eee cee mcm mw me em cc ee wie ieee a 


cee whem mee me eer mr we ieee eee ee 


Neen ee ee ne a nn etiam einem 


Oeol? 356-00 
CeLOOOCE Ol 


O.81737E CO 
Oe-lLOOOCE O1 


0.43045E QO 
0.59193E 00 
0.78334E 00 
Q-l1QOQ00E Cl 


0.99654E-01 0.30318E-0O1 0-62289E-01 


Q0.66789E-O1l 
0.11359E OO 
O.1L7496E 00 


QO-12667E QO 


ccc eee eee ei = 


0.19886E O00 
Q0-28276E CO 


a ee ce ee lee 


ere ee ec ee ee es ae S| a ee 


0-60 0.10 
0.70 0.10 


0.31459E O00 
0.42232E 00 


0.31999E 00 
0-.42507E 90 


0.25489E QO 
0.35664E 00 


0.38026E 00 
O.49115E CO 


0.54782E 00 0.54916E 00 0.48231E 00 0.-61315E CO 


cme wee ct ee 


ee ee ee ee ee ee eee eee eee ee es eee 


‘em cme wh ee eee ee ee ee ee eee 


ec ee ee ee eee eee SS ee ee ee ee es 


0.12058E 00 0.18421E 00 0.86230E—-01 0.16354E CC 
O.23051E 00 0.14186E CO 0.24837E CO 
0.27086E 00 0.25536E 00 0.21022E 00 0.33975E CO 
0.36380E 00 0.37836E CO 0.29428E 00 0.43903E CC 


comm eee 


—— oe ie ae - = 


0.58918E 00 0.5S381E 00 0.51818E 00 0.65874E 


C.19044E 60 


0.47CO7E 00 0.47845E 00 0.39643E 00 0.54594E 60 


—_— —— EE 


= — -—  —— 


00 


O.7L931E 00 O.72i171E 00 0.660COE 00 O.77452E CO 
0.85752E 00 0.85853E 00 0.82120E 00 0.88960E OC 


mee ee eee em ws Ss es es ws - ee ss 


mem ee eee es ae es Ss Se ee sr ere 


-197- 





—_ mec cc ec ww ee ee i ee 


et eee ee eee 





0.30 0.15  —0.21924E 00 0. Re 8673E 00 0. e 16304E 00 Oy 28544E CO 
0.40 0.15 0.30495E 00 0.34939E 00 0.23632E 00 0.38193E CC 
0.50 0.15 0.39994E 00 0.42846E 00 0.32308E 00 0.481S9E CC 
0.60 0.15 0.50473E 00 0.52251E 00 0.42516E 00 0.58551E CO 
0.70 0.15 __ (0 .61887E 00 0.62952E 00 0.54380E 00 0.69132E 00 
0.80 0.15 O.74100E 00 0.74692E 00 0.67948E 00 0.79739E 00 
0.90 0.15 ___ 0. 86897E 00 0.87159E 00 0.83190E 00 0.90120E 00 
1.00 0.15 0. 10000E O1 O.1L000CE 01 O.10000E OL O.10000E O01 
G10 0.17 0.77320E-01 0.26988E 00 0.52922E-01 0.10873E O00 
0.20 0.1? 0.15685E 00 0.2S5S607E 00 0.11180E 00 0.21204E 00 
0.30 0.17  _0.24059EF 00 0.33922F CO _0.17866E 00 0.31279E 0C 
0.40 0.17 0. 33013E 00 0.39853E 00 0.25549E 00 0.41291E 00 
0.50 O17 0442650E 00 0.47286E 00 0.34412E 00 0.51339E 00 
0.60 0.17 0.53006E 00 0.56061E 00 0.44606E 00 0.61429E 00 
0.70 0.17 0.64048E 00 0.65979F 00 0.56236E 00 0.71491E 00 
0.80 0.17 0.75674E 00 0.76798E 00 0.69355E 00 0.81390E CC 
0.90 O17 O487725E 00 0.88240F 00 0.83961E 00 0.90955E CO 
1.00 0.17 O.LOOOOE O1 0.10000E O1 C.LOOOOE O01 O.1O000E O01 
0.10 0.20 0.83324E-O1 0.32281E 00 0.56958E-01 O11 7C2E CO 
Mec 0.20 0.16828E 00 0.34748E 00 0.11980E 00 0.22720E 00 
0.30 0.20 0.25633E 00 0.38803E 00 0.19012E 00 0.33285E 00 
0.40 0.20 0.34866E 00 0.44346CE 00 0.26953E 00 0.43560E CO 
0.50 0.20  0.44600F 00 0.51297F 00 0.35951E 00 0.53634E 60 
0.60 0.20 0.54864E 00 0.59455E 00 0.46121E 00 0.63529E CO 
0670 0.20 __ 0465630E 00 0.68643F 00 0.57587E 00 0.73209E 00 
0.80 0.20 0. 76825E 00 C.78634E 00 0.70378E CO 0.82590E 00 
0.90 0.20 0.88330E 00 0.89176E CO 0.84521E CO 0.91561E G0 
1.00 0.20 O.1LOQ000E O1 0.10000E 01 C.1LOOCOE O1 O.10000E 01 
O10 0622 __ 00 87745E-O1 0.37243E 00 0.59914E-O1 C.12310E CO 
mee0 0.22 0.17669E 00 0.39549E 00 0.12565E 00 0.23831E 00 
0-30 0.22 O.-26791E 00 0.643333EF 00 9-19851E 00 0.34754E 00 
0.40 0.22 0. 36228E 00 064851CE 00 0.27980E CO 0.45220E OC 
0.50 0.22 0.46033E 00 0.5496CE 00 0.37075E 00 0.55311E CO 
feeO0 0.22 0.56227E 00 0662531E 00 0.47245E 00 C.65063E CC 
0.70 0.22 _ 0. 6679CE 00 0.7103SE_ 00 0.58574E_00 0.74463E 00 
0.80 0.22 0.77668E 00 0.8C276E 00 0.71124E 00 0.83466E CO 
0.90 0.22  __0.88774E 00 C.9CCLIE 00 0.84930E 00 0.92004E GO 
1.00 0-22... O-1LOO000E 01 O.LOOGOCE 01 O.1LOOGOE Cl C.10000E O1 


— oe oe 
i eel el 
emcee ea 


ss eee cc cc rc ew i 
eer ee eee ee 
em ee ea 


— — ee 
—_ ces eee ee ee ee ee ee ee ee ee eee ass oo ee we OE Oe” 
ee ee 








CHART l 
x * % x x *x bg Bs K 4 x 
1.000 * 49 * 
0.900 * l 2 3°37 x 
(e000 +*+ #  }.}.®©}®}©.. 1 2 3°45 79 ee * 
e700 * | re NE MO = 
ias00* +t: °° |; 2 3 4 5 6789 " ~ 205 eee 
GO2000 * 1 2 3 4 56 78S x 


ew re mm i i ce we wwe i = = 


eee me ee rere wr aes es 


Sep eere 2 8 34 5 6189 ee eee a 
0.280 *1 2 3 4 56 89 r 
0.180 *123 4 579 ie 
Oo100 *23579 0) A ty 
x x x x xx x x x * x x 


PPOnmo.1C 0.20 0.30 C.40 0.50 0.60 OvromoOlne0 6.90.0 
aioos 





——— ee ee — 


— 


MATRIX 1 10 ROWS 10 CCLUMNS 
STORAGE MODE 0 GROUP » 1 
-— | COLUMN bE EEE DO 
ROW 1 G.100000E 01 0.1000COE 01 0.100000E 01 
ROW 2 0.900000E CO 0-689206E 00 0.777299E CO 
RCW 3 0.800000E OC 0.4236S6E 00 0.571601E 00 
p__ ROW  _% Oe TCOOODE OG _ __—06229394E 00___ _0.396123E 00 _. 
ROW 5 0.600000E 00 0.109590E 00 O0.257871E CO 
~ ROW 6 __ _C.S5COOO0E 00 __ 9.323286E-C1l_ __—«- 041564418 CO _ 
ROW 7 0.400G00E CC C.163907E-01 C.896093E-01 
ROW 8 0.300000E 00 0.505292E-02 0.47477 1E-01 
ROW 9 0.2C0000E 00 0.137341E-02 0.229545E-01 
_--ROW_ 10 __O-LOGOOOE OC _ __ 06 305951 E-03____0.903708E-02 _ 
mee Rix ____i______J pve eee 10 COLUMNS ~____ 
STORAGE MODE 0 GROUP » 2 
8 AGE ES re 2. it Si. ee 
ROW 21  __O-1Q0000E 01 __O.100000EF O1 ___O-100CCCE OL _ 
ROW 892 0.817346E 0C 0.841358E 00 0.857519E CO 
ROW 3 0.644126E 00 0.688840E 00 0.719308E CO 
ROW 4 0.488323E OC 0.547825E 00 0.589176E 00 
ROW 504 35537TE OO _ __— 0 H22325F OC __ _O-470072E 00 _ 
ROW 6 0.247616E OC 0.314594E 00 0.363798E CO 
RCW _ _7___ 06 16462CE OC 0 225015E OC _ _ Ce 27OB8S6E 00 _ 
ROW 8 0.103279E 00 0.152181E OC C.190444E CO 
ROW 9 G.590816E-901 C.931960E-01 0.120577& 00 
ROW 10 C.265891E-01 0.4404S8E-01 0.583032E-01 
“MATRIX 1 10 ROWS 10 COLUMNS ss 
— STORAGE MODE O GROUP +» 3 


rs pgp pp ee SS ene ——ene eee acne ne 


COLUMN 7 8 9 
ROW Lt O.100000E Ol 0.1000COE 01 O.100000E Cl 
ROW 2 _0.868967F OC ____0.877253F 00 __ _0.883305E CO _ 
ROW 3. £.0.741004E 00 0.756742E 00 0.768246E 00 
ROW 4 0.618868E 00 0.640481E OC 0.656300E 00 
ROW 5 0.504727e OC 0.530061E OO 0.548636E CO 
ROW 6 C.399944E 00  _04426457E OC __0.446003E CO 
~ ROW 7 £02e304955E CC 0.330128 OO 0.348656E CO 
ROW 8 0.219237E 00 0-240588E 00  _04256331E CO _ 
mer aw6CUS)St”*é‘i CLL TBE O00 0©6hU0L156850E 00 0.168277E CO 
ROW 10 0.691950E-01 0.773195E-01 QO. 833238E-Cl 
MATRIX 1 10 ROWS _——=—-10 COLUMNS) 


eee ee ee eee eee ee eee eee ee SE ee se eee se ees ee 





——_—_—_—_ em cr wt cr cm cr cr cc yw www ww we ie ee 





COLUMN 10 
ROW 1 C-100000E O01 
meme 6 UCC BBTISGOE OCC ( 
ROW 3 0.77668CE 00 © 
ea Sh CCU GGTIN4SE OC LLL 
RGW 5 0.562272E OO 
ROW 6 0.460333E CC 
RCW 7 0.362277E OO 
Meeees 68 0 267913E CO. 
ROW 9 G.176689E 00 
ROW 10 0.877451E-01 


ec te, a se ee le ee 


coerce rtm mr cc rr cr rc cc rec cc cmc cc wi ec wc ce ee ew eee 


i i ee ae el 
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ee a es GS ee ie SS ee ee ee ee ee ee ee ee i ee ee ee ee ee ee ee ee en ee en en ee en SS ss 


ee, Ee nn in a a a a a > a a eee ew ee ee eee ee ee eee ee ee eee en nes a= =e eae SS |_| CUD 


Som) Sse ee ee ee ee nn a ie a a i oe a a a oe ae ae ee ee ee ee ee ee ee 


i er 
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—_—> ce ce cee ee eee eee eee eee 
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ee 
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ete 





ce ee ee eee ee 


ne - CGS eee a 


MATRIX 2 10 ROWS 10 COLUMNS 
STCRAGE MODE O GROUP , 1 
ar MOmUMND Ml 0 
ROW 1  O.100000E 01 #0.999999E 00 0.100000E 01 
ROW 2 0-SOCGO00E 00 0.689158E OC O.777291E 00 
ROW 3 0.800000E 00 0.423778E 00 0.571596E 00 
_ ROW 4 ~_0.700000EF OC 0.229871E 00  0.396115E 00 _ 
ROW 5 0.600000E 00 0.105476E OO 0.257800E 00 
___ROW_ 6 __O0.5C00COF 0G 0.456947E-01 —_0.157242E 00 
ROW 7 C.4CO00G6E 00 0.166054E-01 0.897810E-O1 
ROW 8 0-300Q0000E OC 0.498009E-02 0.480379F-O1 
ROW 9 0.200000E 00 0.356913E-02 0.243909E-Cl 
__ ROW 10  _O.100000E 0C 0.373244E-03 001 27263E-01 _ 
_MATRIX == = 2 10 ROWS LO COLUMNS 
STORAGE MODE O GROUP , 2 
— Cl PoeMN 2224 oS CU eee 
___ROW i  __O-iO0S000E Ol  O0-100000E 01  O.100000E O1 | 
ROW 2 C.817368E 00 0.841594E 00 0.858528E 00 
ROW 3 0.644190E 00 0.689454E 00 0-721 7C5E 00 
ROW 4 0.488495 OC 0.549163E 00 0.593811E CO 
___ROW 5 02 355814E 00 = =— 0 425069F O00 «06 4 78449E 00 | 
ROW 6 0.248729E 00 0.319994E OG 0.378359E CO 
___ROW 7 _—6.367068E 00 0.225234F 0G 0.295359E 60 _ 
ROW 8 O0.108085E OC C.170740E 00 0.230511E CC 
ROW 9 0.702469E-Cl 0.125980E 00 0.184211E CO 
ROW 10 0.487320E-01 0.996541E-01 C.156484E 00 
Mmayeix 2 #42110 ROWS 10 COLUMNS ©. 
meee). 6©60OS—téi‘<i‘ SS TQRRAGE MODE OF °§}©GROUP,30—lmté<“‘ 2 2W®”*~=‘i‘S™SS 
COLUMN 7 8 9 
ROW 1 £O.1L00000E 01 #4£44G.1CCCOO0E 01 °&#&«&+0.100000E 61 
> ROW 2 £#0.871591E 00  0.882398F OC 0.891763E 00 _ 
ROW 3 0.746924E 00 0.767981E 00 0.786335E CO 
ROW bs 0.629523E 00 0.6597S4E 00 0.686427E CO 
ROW 5 0.522506E OC 0.560613E 00 C.594554E 00 
RCW 6 0.428458E OC 0.472858E 00 0.512968E 00 
ROW 7 0.349387E 00 0.398534E OC 0.443600E 00 
ROW 8 0.286730E 0G _ 0.339217E OC — 0.388030E 00 _ 
ROW 9  0.241436E CO 0.296071E 00 0.347479E CO 
ROW 10 O-214071E CO C.269883E 00 0.32280SE 00 
MATRIX 2 10 ROWS 10 COCURNS 3 ee 
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mcm ees se ee ae — a — 








STORAGE MODE O GROUP , 4 
COLUMN VG ean ee 
ROW 1 O.100000E 01 
ROW 2 _(0.900113E 00. 
ROW 3 GBesc2758t OG ## §| Bette te 00a ae 
moe = CCH LO0SB89E COU 
RCw 5 0.625308E OC 7 
ROW 6 0.549605E 00 
ROW 7 0.485104E OC 
Meee OO 495S431E OC 
ROW 9 0.395486E CC 
RCW 10 0.372433E 00 
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MATRIX 3 10 ROWS 10 COLUMNS 
STCRAGE MODE O GROUP , 1 
a Gen # tf Sy 0 == Sa oe 
ROW 2 O.1LOO0000E 01 #£40.100000E 01 #0<1000COE Ol 
ROW 2 0.900000E 00 0.661296E 06 0.745336E CO 
ROW 3 0-800000E OCG 0.390233E 00 0.525795E OO 
___ROW 4  _—_O0.700000E 00 0.202861E 00  0.349673E CO _ 
ROW 5 0O.600000E 00 0.930687E-G61 0.218502E CO 
___ROW 6  _0.500000EF 00 0.264024F-01 — 0.127269E CO | 
ROW 7 0.400000E OC 0.128419E-O1 0.699988E-Cl 
ROW 3 0.300000E GO C.380260E-02 0. 356159E-O01 
ROW 9 0e2CO0COE O00 0.993013E-03 0.165389E-Ol 
___ROW_ 10  —0-100000F 00 0.212133E-03 0.625416E-02 | 
memetRIX 3. __ _10 ROWS  _ 10 COLUMNS 
STCRAGE MODE O GROUP , 2 
— SC) | i i Seen 
___ROW 1 = _O.1COO000E Ol  O.100000F O01  0.100CCOE O01 
ROW 2 0. 783344E OC 0.806019E 00 0.821203E OO 
ROW 3 O.59L927E AC 0.632494E 00 0.65SS98E CO 
ROW 4 Oa4seeave OO 0.482313E 00 0.518i77E O00 
_ ROW 5 = —0-300572E 00 0.4356637F OC 0.396426E 00 
ROW 6 0.200998E C 0.254890E OC 0.294283E GO 
___REW 7  O.128274EF 00 O.174963E OC 0.210224E 00 | 
ROW 8 0. 772659E-01 0.113586E CC 0.141865E CO 
ROW 9 0. 424448E-01 0.667886E-O01 0.862302E-Cl 
ROW 10 C.183465F-O61 0.3031 76E-Ol 0.40041 3E-Cl 
"MATRIX 2 #210 ROWS | 10 COLUMNS .}.}.}}©. 
—  §FORAGE MODE 0 | GROUP » 3000 )>) >) 0a 
COLUMN 7 8 9 
"ROW 1 Ocl000GO0E 01. £O-LOOQ000E 01 °#£0.100000E C1 
__ ROW 2  04831905E 00 0.839613F 00 ___0-.845214E 00 _ 
ROW 3 0.679485E OC 0.693550E 00 0.70378CE CO 
ROW 4 0.543800E 00 0.562358E OC O.575874E 00 
ROW 5 0.425159E 00 0.446060E 00 0.461308E OO 
___ ROW 6 O0.323077E 00 0.344124EF 00  _0.3595CSE_ CO_ 
ROW 7 £«.66236322E 00 0.255493E OC 0.269533E OO 
ROW 8 0.163038E OC 0.178660E 00 €.190122E 00 _ 
ROK 9  °&£«0.100923e 00 O.1LLISO4E OO O.1L19797E 00 
ROW 10 0.474349F-01 0.529225E-Cl 0.569578E-01 
MATRIX 3 i0 ROWS 10 =C OU UMN 
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COLUMN 10 
ROW 1 C.1COO0GOE O01 
Meee eC UO ONOZGBE CCC 
ROW 3 Os721242£ 00 }°.}©=©)>L.c 0 5 5 nen 
moe 4 CCC SB5739E 00 rt—~—C‘CSC‘C*S 
ROW 5 0.472445E 00 (400 
ROW 6 C.370754E CO 
ROW ? 0.279803E 00 
Meee = O61 DSSI2ZE OO ee 
ROW 9 Q.125652E 00 7 — 
ROW 10 C.599143E-01 
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x x % ad * aK an x x x *x 
1.00C * 69% 
— rr 
— “22. 
oe. enene 
a a 
0.500 * 1 2.3 __¢*_)5) $6137 eee 2 
Meee) 2 38 4 8 IBS ee =f 
Ss oe 


jm wee ee es es es es es es ee ces es ee es es es ee ee es es es ees es ew cs i es ee ieee ees es eee en ee ees ees eee Re eee 


emp me ee ew es ee we ee sn nn eae es ae 


mee seme wi mc ee ee ee eee 


mem lem cme me ee ee ee ee ee eee ea ss a 8 Oe eS ee _ ss se 
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MATRIX 4 10 ROWS 10 COLUMNS _ 
pe ee ee ee 


STORAGE MODE O GROUP , 1 

— | COUUMN OQ SB 
ROW 1 O.100000E 01 0. 100000E Ol C.100000E O1 

ROW 2 0.900000E OC 0.716376E 00 0.8074L6E 00 

ROW 3 0.800000E 00 0.457947E OC 0.617032E CO 
___ROW 4 ~_0-?700000E OC  _0.257891E 00 04 444528E 00 _ 

ROW 5 0.600000E 00 O.128171E O00 0.300911E 00 
___ROW_ 6 == 0.500000E CC  0.393896E-01 —_0.189869E 00 _ 

RCW 7 0.400000E 00 0.207545E-01 O.113127E CO 

ROW 8 0.300000E 00 0.665832E-02 0.623550E-01 

ROW 9 GO. 200000E OC 0.188375E-02 0.313673€-O1 
___ ROW 10  _©.1Q00000E 00 0.436425E-03  __0.128500E-01 _ 
MATRIX = _ 4 10 ROWS dO COLUMNS) 

STORAGE MODE O GROUP ,» 2 

___ a BUCURNe 4 eee 
__-_ROW__ 1 = 0.100000E O01  O.-100000E 01  C.100000E 01 | 

ROW 2 0.848590E OC 0.873154E 90 0.889603E CO 

ROW 3 0.694639E 00 0.742246E 00 0.774522E 00 

ROW 4 0.547213E 00 C.613148E 00 0.658741E 00 
Moen = A 413934E 00 =—0 491145E 00  §=—0.545941E 00 

ROW 6 0.299862E OC 0.380262E 00 0.439030E O00 
ROW __7 ___0.207307E OC___0.282761F OC____0.339748E 00 

ROW 8 0.135273E OC 0.198860E 00 0.248368E 00 

ROW 9 0.804995E-021 0.126669E OC 0.163541E 00 

ROW 10 0.376939E-O01 0.622888E-01 0.822664E-O1 
“MATRIX 4 4120 ROWS | 10 COLUMNS ©. 
meen ))6|)™SCUSTGRAGE MODE O©=F”)™D™™—CG REP 5 D~~™t~é<“‘OSOSCSOOCOCO:C*~*”W 

COLUMN 7 8 9 

— ROW 1 O.10000CE 01 £O-.100000E 01  °}#&£2©0.100000E O01 
Mowe 6=—CU PC COOL IOGE 00 === GOS D4S6E OO  __O.9LS6T3SE 00 

ROW 3 0.7S7390E 00 0.813896E O00 0.825901E O00 

ROW 4 O.691315E 00 0.7149C8E 00 0.732090E OO 

ROW 5 0.585510E OO 0.614294E 00 0.635293E CO 

ROW 6 0.481987E 00 0.513387E 00 0.536338E OO 

ROW 7 0.381926E CO 0.412908E 00 0.435598E 60 

ROW 8 0.285437E OC  __0.312787E 00 __0.332854E CO 

ROW 9 °&£06191407E OO 0.212042E OC 0.227202E 00 
ROW 10 0.974566E-01 0. 108731E 00 | Gomero77 econ 
MATRIX 4 10 ROWS 10 COLUMNS si 
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COLUMN 10 

ROW 1 C.100000E O1 
_______ ROW 2 = 0.920038E OC 

ROW 3 0.834659E 00 © 
eee = 0 TS4O31E 000 =e 

ROW 5 0.650631E 06  ...}.}.}}}©= ©" 

ROW 6 eS53115E OC 

RCW OO? 0.452196E 00 

ROW 8 0.347543E 00 


—_—_—_e one mc we rr wr rr crm cr rw mc me rc ome ee es es es eo i 


ROW 2 O0.238306E OC 
0 0.123096E 00 
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APPENDIX A5 


SUBROUTINES USEFUL FOR COMPUTING LPCM STEADY-STATE SOLUTION 

For completeness and for the convenience of the reader, this 
appendix presents a listing of the Fortran IV statements in the sub- 
routines useful for computing the LPCM steady-state response. As 
mentioned previously it is expected that there are easier and more 
efficient methods which could be developed and used to find the 
steady-state for special cases of the LPCM, These subroutines will 
solve the general LPCM steady-state with the subroutine AFCT given 
here, Several of the subroutines here are also used directly in PLPBV, 
namely EIGEN and GELG. 

Four of the subroutines used here are taken, less comment cards, 
directly from reference: (I-3), the Scientific Subroutine Package. 
These are subroutines LBVP, EIGEN, GELG, and LOC. The remaining 
three subroutines, AFCT, FCT, DFCT, were written specifically for 
the LPCM. One additional subroutine, OUTP, is required for operation 
of the subroutine LBVP and must be furnished by the user. All of 


these subroutines compile under WATFOR in a total time under 20 seconds. 
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CSUBROUTINE ULBVP(PRMT RC Ry Yy,NERVeNNDIM, THLE, AFCT,FCT,DFCT, WTP, 
TAUIX,A) 

DEMENSTON PRMT (1), BOLI eCO1»RCL I »VE1) sDERY(L) pAIIX( 24,1) ,A01) 

TF CPRMT (39% (PRMT(2)-PRMT(1)))7,1,3 

} YTHLF=12 
RFTURN 

BetHLF=13 
RFTURN 

3 KK=-NNIM 
I B=0 
TCc=0 
DQ 7 K=1,NDIM 
AUX(15,K) =DERY(K) 
AUX(1L,K)=1. 
AUXC1L7,K)=1 e 
KK=KK4+ND IM 
DP 4&4 T=1,sNDIM 
TYT=KK+I 
eC S0CT1)) 5, 4,5 

4 CONTINUE 
TR=TR+1 
AUX(1L,»K)=H°. 

5 ON 6 T=1,NDIM 
TT=KK+I 
TF(CCTIIV)7,6,7 

6 CONTINUE 
fei C+ 
AUX(17,K)=C6. 

7 CONTINUE 
me. to- TR j}Re 11,11 

RB H=PRMT(2) 
PPMT (7? })=PRMT(1) 
PRMT(LI=H 
PRMT(32}=—-PRMT (2) 
DP 9 T=1,NDIM 

9 AUXCI7,TI=AUX(1,1) 
TI=NOIM*NDEM 
DN 19 T=1,IT 
H=R(T) 

B(T)=C(T) 
1% C(T)=H 
TY X=PRMT (2) 
CALL FCTIX,Y) 
CALL NFCT(X, NERY) 
NA 12 I=1,NOIM 
AUX(18,TI=V(T) 
12 AUX(19,S)=DERY(T) 
K=0 
KK=0 
100 K=K+} 
TFCAUXC 17 KIITCS,198,191 
101 X=PRMT(2) 
CALL AFCT(X,A) 
SUM=0 < 
GL=AUX(18,K) 
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102 
193 
194 


FOS 
104 


Le 7 
108 


LO9 
a 


el 2 


PL 3 


114 


Li6 
il 7 
aon 


PO! 


2C? 
2X3 


DGL=AUX(19,K) 

TT=K 

DO 104 T=l1,NDIM 
H=-A( TT) 

DERY(T)=H 

AUX 26,7 )=R(T) 
Y(T)=e, 

ITF( 1-K) 193,102,103 
Y(T)=1. 
AGL=EDGL+#H*AUX(1LA,T) 
TT=FT+NDIM 
XFND=PRMT (1) 
H=.%625*(XEND-X) 

IT SW=9 

GOTN 409 

TEC THLF-1%) 106,106,117 
AD LOT T=1,NDIM 
KK=KK+ 1] 

H=C(KK) 

R(7T)=AUX (20,1) +H*SUM 
Tr=! 

DN 107 J=1,NDIM 
R(TTI=B( TT) +H*Y( JS) 
TY=TTtnprM 

GOTTA 1N9 

KK=KK+NN TM 
IF(K-NDIM)100,110,110 
X=PRMT(4) 

CALL GELG(R,ByNDIM,1,X,7T) 
PECL it sll 26132 
THLFH=14% 

RETURN 

PRMT(5 _—o. 

IHLF=-T 

X=PRMT (1) 
XFND=PRMT (2) 

H=PRMT (23) 

DM 113 T=l,NDIM 
Y(T)=R(T) 

TSw=] 

I SW2=12 

GOTO 2AM 

I SW3=-1 

GOTTA 390Nn 
JFCTHLF)490,4900,117 
RF TURN 

CALL AFCT(X,A) 

IF (CT SW)201,701 92705 


LL=6 

NN 203 M=lL,NDIM 
HS=0. 

09 29? L=1,NDI™M 
LL=LL+1 
HS=HS-A(LL) *Y(L) 
DERYCMY=HS 


2m 3 





204 GOTO(502, 504450644074 41594185608 9617 963296345421 9115) 9 1 SW? 
mee CALL FCT(X,NDERY) 
DQ 207 M=1,NDIM 
LL=M-NDIM 
Mo= 0's 
DY 20h L=l,NDIM 
LL=LL+NOEM 
206 HS=HS+A(LL) *VOL) 
2OT DERY(MJ=HS4+DERY (M) 
GATA 204 
ANC TE( TSW) 301,301,395 
30) CALL FCT(X,R) 
GU=9, 
NGU=90. 
NN 302? L=l,NDIM 
GU=GU4+Y(L)*RQ(L) 
399 DGU=DGUFDERY(LI*RI(L ) 
CALL DFCT(X,R) 
DA 303 L=1,NDIM 
393 DGU=ENGUEV(IL IR (LD) 
SUM=SUM+ .5*H*((GL4+GU) +. 1566667*H*( DGL-DGU) ) 
GL=GU 
OGL=DGU 
394 T° (1TSW3)116,422?,61 8 
395 CALL OVUTP(X,Y,DERY,FHLFESNDIM,PRMT) 
TF(PRMT(5)) 117,394,117 
409 N=] 
XST=X 
THLF=9 
DN 401 T=H=1,NDTM 
AUXC1LA,F)=C. 
AUX(1L,TI=YCT) 
407 AUX(8,T)=DERYC(T) 
TSwWi =] 
GOTO 4590 
402 X=X+H 
N99 493 T=1,NNIM 
493 AUX(2,T)=¥CT) 
404 THLF=THLF+] 
X=X-H 
DM 405 fT=1,NDIM 
405 AUX(4,F)=AUX(2,7T ) 
H=.5%*H 
N= 1 
1SWl=? 
GNTQ 50 
406 X=X+H 
ISW2=4 
GNTA 200 
4QN7T N=? 
DM 408 T=1,NDIM 
AWIX(2, TY=YCT) 
4098 AUX(9,TI=NERY(T) 
TSW1=3 
GNTA 890 
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499 DON 414 T=1,NDIM 
7J=ABS(YCT)) 
TFCZ-1.9410 5411541} 
410 7Z=1. 
41] DFELT=.96666667* ABS(Y(T)-AUX(4,1T)) 
IF (TSW) 413,413,412 
412 DELT=AUX(15,1)*NEI T 
413 TF(DFL T-Z*PRMT (4) 414,414,479 
414 CONTINUE 
X=X+H 
TSW2=5 
GOTTA 200 
415 DN 41464 JT=1,NDIM 
AUIX(3 ,TI=Y(T) 
416 AUX(1C,T)=DERY(T) 
N=3 
ISW1=4 
GNTN 5a4 
417 N=] 
X= X+H 
1SW2=4 
Gatm 200 
418 X=XST 
NN 419 T=L,yNNITM 
AUX(11,TIY=DERYCT) 
AZLOCYCTP=AAUK() ~TPF#H*OE 3 75*AUX (BRT) +. TIOLE667TXAUX(9, T) 
1- «20832333 *AUX( 10,794.94) 66667*XDERY(T) ) 
42) X=X+H 
N=N+ ] 
TSw2=1} 
GOAT 2700. 
421 TSwr=f 
GOTO 300 
422 TF(N—-4)423,609,69%C 
423 0N 424 Y=l,NDIM 
ANIX(N,TI=HV(T) 
424 AUX(N47, TI=DERYCL) 
IF (N-329425,427,600 
425 09 424 Y=1,NDIM 
DEL T=AUX(9, T) #AUX(9,1T) 
DELTHNEL T+DELT 
4926 YC TI=HAAUKO1, T) +. 3333333*H(AUX(BgTI+DELTHAUX( LC, T)) 
GOT 42" 
427 DN 429 T=1,NDIM 
NDFEL T=AUK(9,T)#AUX(19,T) 
DFLT=NEt T4#DFLT+DELT 
4P7R YC TI=AUX(1 LT) to 2 75*#H*(AUX( 8B, TP) FNEL T+HAUX(1T 1, TD) 
GNTN 42C 
429 YTECTHEF-19) 404,430 ,43% 
4320 THLF=11 
X=X+H 
TE(TSW)195,1905,114 
500 P= xX 
NN 5901 T=1,NDIM 


X=H*X*AUX(N+7, 1) 
Se 





AUX(5,1T)=X 
S°r YC TI=AUXIN, 194.4 *%X 
X=Z7+.4*H 
T SW2=] 
GNTA 20¢ 
5®2? DNA 503 T=1,NDIM 
X=H*DERY(T) 
AUX(6,1)=X 
5V3 VC TI=AUX(N GT) 4.2969 776XAUXK(5,1)+.1587596*%X 
X=74+.45573727%H 
[SW2=2 
GNTN 200 
594 DOM 595 T=1,NDIM 
X=H*NERV(T) 
SVS VCT)=AUX(N, T) 4.2L RLID4*AUX (5, 19-3 .NSDIGSEAUX (651) 43.8 32865*X 
X=Z7+H 
T SW? =2 
GNATO 2PO 
506 NN 507 T=1,NDIM 
SOUTOY(TP=AUX(N, 1) 4+.1747603*AUX (5, 19-.5514 B07 XAUX (6,7) 
Ltl.2905536¥%AUX(7, 1 )+#.1 711 848*H*XDERY( 1) 
X=7 
GNTN(402 ,406,409,417), ISwW! 
HNN TSTFP=3 
AOY TFIN-8)604,607 ,604 
62 DN 6N3 N=2,7 
DA 693 T=1i,NDIM 
AUX (N-1,T)=AUX(N,T) 
60% AUX(N+6,1T)=AUXCN+7, 1) 
N=7 
604 N=N+1 
AN 605 T=1,NDIM 
AUX(N-1,7T)=Y(T) 
605 AUXIN+6, 1J=DERY(T) 
X= X+H | 
696 ISTFP=ISTEP+H1 
NN 607 T=1,NOIM 
ODELT=AUX(N—4, 17 341.3333 33 *H* (AUX (N+6,1)+AUX(N46,1T)-AUX(Nt5, I) + 
TAUXK(N4#4, 1) +AUX(N44,7)) 
YC FI=DNEL T-. 9256198*AUX (16,7) 
6C7 AUX(16,1T)=DFELT 
I SW? =7 
GOTA 200 
698 DON 699 T=1,NDIM 
ONFEL T=.125%09. AUX (N—1 9 I) AUX (N-3, 779 +3 .eH*E(DERVOT) FAUX (N46,7) + 
AUX(16,1)=AUX(16,7)-DELT 
6099 YCL) =NELT+.9 743803 7¥AUX(16,T) 
DEL T=, 
NO 616 T=1,ND7M 
Z=ARSCY(T)) 
TF(Z-1.3610,611,611 
610 7=1. 
611 Z=ARS(AUK(16,1)972 
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TE(ISWI613,613,612 
617 Z=AUX(15,1)*Z 
613 TF(Z-PRMT(4))614,614,628 
614 FFC(NEL T-Z)615,616,616 
615 DFLT=7 
616 CONTINUEF 
TSW2=8 
GNTOQ ?70C 
617 TSW3=1 
GOTO 30° 
618 FECH*(X—XENND) 9619 ,621,6271 
619 TFCABS(X~XEND)—.1*ABS(H) 1671 ,629,620 
62) TF(DELT—-.N2*PRMT (4) ) 622,627, 6} 
671 TFC(ISW)195,105,117 
Ore? 1 F( 1THLF)601,601,623 
673 TFE(N-7)601,674,674 
674 TEC ISTEP-4)4601,625,6275 
4625 IMNN=ISTEP/S2 
IFC TSTEP—IMON-IMND)601 ,6276,69) 
626 H=H+H 
THLF=THLF-1 
ISTFP=0 
DN 627 T=1L,yNDIM 
AUX(N=-1,T)=AUX(N-?7, 1T) 
AUX (N-2,1) =AUX(N-4, 1) 
AUX (N—32,T)=AUX(N-6,T) 
AUX (N4+6,1)=AUX(N45,17) 
AlIIX(N4+5,1)=AUX (N43, 7) 
AUX(N+4, 10)=AUX( N41, 1) 
NEL T=AUX (N46 ,1) + AUK ON45, T) 
DELT=DELT*+DELT#ANELT 
62 TOAUK(16,T) =8.9629 62*(VY (1) -AUX(N-3,7) 9-3-3261 1L1L1*H* (DFQAYCT)4DEL T 
1+AUX(N+4,7)) 
GNT 491 
628 YHLF=ITHLF+1 
TEC THLF-10)630 463% ,6279 
629 TF(1SWI105,105,114 
6290 H=,.5*H 
ISTEP=0 
DO 631 T=1,NDI™M 
OY (T)=.00390625*(80. *AUX(N—1_9 1) +135. *AUX(N-2, 1) 440. XAUX (N-3,1) + 
TAUX(N—4,T))-—- 11718 75% (AUX (N+t6,7)-6.*%AUX( N45, TI-AUX(N4+4,1)) *H 
QAUX(N-4,7T)=2-90390625*(17.*% AUX(N-1, 1) #1325. *AUX (N-2,17)+ 
11°08. *AUX(N-3, TD) + AUXIN—4,_) 1) )-207343 758 (AUX (N+6,1T) +18. FAUX (N45, 1)- 
29.* AUX (N4+4,70))%*H 
AUX(N-3,7)=AUX(N-?,7T) 
S31 AYX(N+4, TF) =ANXK(N45,7) 
NELT=X—H 
X=DEL T-(H+H) 
TSW2=9 
GA™| ?7o0r 
632? DON 633 T=l1,NNDIM 
AUX(N-?2,T3=YCT) 
AUX(N+5, TI=DERYCT } 


632% Y(T)I=AUX(N-4,7T) 
-217- 





X= X- (H+H) 
TSW2=19 
er wot 
34 X=DELT 
NN 635 T=1,NDIM 
NELT=AUXON45,1)4+AUX( N44, 1) 
NELTHNELT+ENELT+NEcy 
CAUX(16,7T)=8-967963%( AUX(N-L,PIKVOT)I—3 361 LILI *XH*( ALIX (N46, TP #NELT 
L+NDERY(1)) 
35 AUX(N4+3,IF)=DERY(T) 
GOTO 696 
END 
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SUBROUTINE GELG(R yA yMyN,FEPS,TER) 
NIMENSTON A(1),R (1) 
TF (M) 23,23,1 
TFR=0 

PTv=N, 

MM=MxM 

NM=NXM 

DN 3 {=1,™M 
TR=ABS( ACL) ) 

mit fh -P(V)3,392 
PYV=TR 

I=L 

CONTINUE 
TOL=EPS*PTV 
LST=1 

Nn 17 K=1,M 

IF (PIV) 23,23,4 
PEL TER) 7,5,7 

IF C(PIV-TOL) 6,6,7 
TFR=K-]1 

PYTVI=1./ ACT) 

J=( 1-1) /M 

[= f-j*x4—K 

J= J+1-K 

DD 8 L=KyNM,M 
LL=L+T 
TR=EPIVIXR(LL) 
R(LLI=HRIL) 
R(L)=TA 
ITF(K—-MI9,18,18 
LEND=LST+M—-K 

PREC II 12512, 1 
TJ=J*M 

NAO 11 L=LST,LEND 
TR=A(L) 

LL=L+Tf 
A(L)=A(LL) 
A(LL)=TB 

DY 13 L=LSTyMMy™ 
LL=L+]7 
TB=PITVI*ACTLL) 
A(LLI=A(L) 
A(L)=T8 

ATLST)=J 

PITV=0. 

LST=LST+] 

J=O 

DN 16 [I=tLST,LFEND 
PIVYI=-A( TT) 
TST=FT+M 

J=JS+1 

PA 15 L=ISTyMM,M 
LL=L-J 
A(LI=A(LI+PIVIXACLL) 


TR=ARS(A(L)) 
-219- 








14 


15 


16 


18 
uo 


21 


223 


IF(TB-PIV) 15,154,114 
PTV=TR 

T=t 

CONTINUE 

DN 16 L=K,NM,M 
LI=L+J 
R(LLI=HRILLIFPIVIXRIL) 
tST=LST+M 

TF OM-1)23,22,19 
IS T=MM#M 

LST=M+4]1 

DO 21 JY=?_M 

f= ST—] 
IST=I1ST-LST 
L=IST-—M 
E=AlI}+.5 

DD AL J=HTT, NMy™ 
Wma=R( 3) 

LL=J 

DN 20 K=IST»MM,™ 
LL=Ltt+l 
TR=TB-A(K)*R(LL) 
K= J+. 

R(JI=HR(K) 
R(K)J=TA 

RETURN 

IFER=-1 

RETURN 

END 
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SUBROUTINE AFCT(X,A) 
THIS SUBROUTINE IS A PORTION OF PLPRV 


THESE DIMENSTONS MUST BF CHANGED IF THE INPUT 
POLYNOMTAL DIMENSTONS OF PLPBV ARE CHANGEN. 
DIMENSTON P (3) ,Q(3) 


MATRIX A MUST BE IN VECTOR FORM 
NOYIMENSTON ACY) 


Py>Q, AND N ARE POLYNOMTAL COEFFICIENTS AND MUST SF 
MAIN PROGRAM AND PLACED IN COMMON AREA S . 


CIMMON /S/P,0,N 
FPSTLON IS ZERO VALUE FOR THIS SUBROUTINE 
FPS=1 .OF-97 
PY=P(1) 
P3=Q(1) 
P4=0,0 
TF CN-1)10,11,12 
PyY=P14+P(7)*X 
P4=Q(?) 
GN TO 19 
PT=P]14+P(?7)*X 
P3=P 3247 ,90*P(2)40( 2 )*X 
P4=Q(2) 
XIMM=1 .9 
AN 20 T=2,_N 
XIM=X J MM*X 
XJ=XIM*X 
RT=FLOATCT) 
PY=PL4EP(T4+1)*XT 
P3=P34? ,.O*XRIXP(T +1) *XIM+O(C T4+1) *XT 
P4=P44RT*#((RT-1.0)*%P( T4l )*XIMM4Q(T 41) XX9M) 
X TMM= XT MM* X 
CONTINUE 
TF(ABS(P1)—EPS)13513914 
A ( 2)=C rae) 
Meo v=o. C 
A(4)=°0,9 
TF(ABS(P3)-—£PS)15,15,16 
A(1)=9.0 
RFTURN 
A(1L)=-P4/P3 
RETURN 
ACV) =F." 
A(2)=-P4/P1 
A(3)=1.f 
A(4)=-93/P] 
RFE TURN 


END 
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SUBROUTINE FCT{X,F) 

THIS SUBROUTINE TS A PORTION NF PLPRV 
DIMENSTON F (2) 

F(L)=°,0 

F(2)=C.0 

RETURN 

END 


SUBROUTINE DFCT(X,DF) 

THIS SURROUTINE ITS A PORTION OF PLPRYV 
DIMENSTON DFO?) 

DEtL)=9.0 

Pee) =o 6 9 

RFE TURN 

END 


SUBROUTINE LOCC T,J,IRyNeMeMS) 
{ X=] 

JX =J 

TECMS-1t) 10,20,30 
TRX=NX( JX—-1 +7 X 

GMO TQ 36 

TFC IX-JSX) 224924424 
TRX=TX+0IX*IX-JX)/? 
GON TO 36 

TRX=IX+C IX*TK-1TX)/2 
GN TO 36 

TR XK=0 

IFC IX-JX) 36532,36 
IRX=[¥ 

IT R=T7 RX 

RETURN 

END 
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SUBROUTINE FCT{X,F) 

THIS SUBROUTINE TS A PORTION NF PLPRV 
DYIMENSTON F(2) 

F(L)=9.0° 

F(2)=C.0 

RETURN 

END 


SUBROUTINE DFECT(X,DF) 

THIS SURRNUTINE ITS A PORTTON OF PLPRYV 
DIMENSTON DFO?) 

WET )= 9.20 

RETURN 

END 


SUBROUTINE LOC CT,J5, TR NyMyMS) 
[X=] 

JX =J 

IF (MS-1) 10,20, 30 
IRX=N*X(SX—-1 +7 X 

GN TQ 36 

TFC IX-IX) 22924924 
TRX=TX+( ISX*IK-UX)/? 
GN TO 36 
TRX=SK+( IX TX-1X)/2 
GM TON 36 

TRX=0 

IFC IX-JX) 36,32,36 
IRX=!¥X 

TR=TR¥X 

RETURN 

END 
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bas. 
an 
25 


30 
35 


4a 


45 
55 
60 


6? 
65 


68 


cad 
75 


78 


gn 
85 


he. 
95 
pAA 


SUBROUTINE EIGEN(CA,R aN MV) 
DIMENSTON AC1L),REI1) 
RANGE=1.0E-6 

IFCMV-1) 10,25,10 

ITOQ=-N 

NN 20 J=ly_N 

[Q=TQ+4+N 

DO 290 T=1,N 

TJ=IQ+] 

R(IJ)=0.0 

teet—s)} 20,15,2° 
R(TJI=1.0 

CONTINUE 

ANORM=9.0 

N97 35 T=1,N 

NN 35 J=I,N 

TF(I-J) 309, 35, 3 
TA=14+(J*J-J)/2 
ANORM=ANORM4A(TA)*A( TA) 
CONT TINUE 

ITF CANAIRM) 165,165, 40 
ANORM=1.414*SQRT ( ANORM) 
ANRMX=ANORM*RANGE/FLOAT(N) 
INN=0 

THR=ANDIRM 

THR=THR /FLOATI(N) 

L=1 

M=L¢] 

MQ=(M*M=M) /2 
LO=(L*¥L-L)/2 

tL M=L_+MQ 

TFC ARSCAC(LM) J-THR) 139,65,65 
PAp=1 

LL=L+L9 

MM=M+4+MQ) 
X¥=9,5*(ACLLI-~A(MM)) 
Y=-A(LM)/S SORT CACLM) XA (LM) +X*X) 
IF(X) 79,75,75 

Y=~-Y 

SINX=V/ SORT(?2.0*(1.04( SQRT(1.90-Y*Y)))) 
SINX2=SINX*®SINX 

COSX= SQRTCL.QO-SINX2)} 
CNSXK 2=CNSX*COSX 

STINTS =SINXK*COSX 
TLQ=N*(L=1) 

IMQ=N* (M-1) 

D9 125 T=1,_N 
1O=(1*{-1)/? 

Te Ci—t ) 80,115,809 

IFC I—M) 85,115,908 
IM=Y+MQ 

GN TA 95 

1M=M4+1Q 

IFC T=-~L) 190,195,195 


TL=I4+LQ 
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GO Tn 11% 

IL=L+IQ 
X=A(ITL)*COSX—AC TM) *S INX 
ACIM)=AC IL) *SINX#A( TM) *COSX 
At It) =x 

TFECMV-1) 120,125,120 
TLR=ILQ+E 

T MR=IMQ+T 
X=R(ILR)*COSX-RE IMR) *SINX 
R(IMRI=HRIILRIXSINX4R (CT MR) *COSX 
R{(ILR)=X 

CONTINUE 

X=2.0*A( LM) *SINCS 

Y=A(L1I )*CNSX24+A( MM) XSTNX?-X 
X=AC(LLI*SINX24+ AMM) *COSX 74X 
AC LM) =CACLLI—ACMM)) XSTNCS#A( LM) *(COSK2-SINX2) 
A(LL)=Y 

A( MM) =X 

ITF(M-N) 135,140,135 

M= M+] 

GN Ta 60 

IF(L—(N-1)) 145,159,145 
L=L+1 

GN TON 55 

TFC IND-1) 160,155,160 

IND=0 

GO 10 50 

TFC THR—-ANRMX) 165,165,495 
IOQO=-N 

DA 1ABS Y=1,N 

TQ=I1Q+4N 

LL=I+(T*{-1)/2 

JQ=N%( 1-2) 

DN 185 J=I,N 

JQ=IQ+N 

MM=J+(J*J-4)/2 
IFCACLLI-A(MM)) 170,195,185 
X=A(LL) 

A(LLY=A(MM) 

A(MM)=X 

TF(MV-1) 175,185,175 

DQ 189 K=1,N 

TLR=TO+K 

TMR=JQO+K 

X=R( TLR) 

RCFILRI=HR( TMR) 

R{ IMR) =X 

CANTINUE 

RETURN 

END 
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APPENDIX A6 


SUBROUTINic PLPBV 

This appendix presents a complete listing of a suggested format 
and some of the computation statements necessary to solve the PLPBV 
model which is a special case of the LPCM, As stated in the subroutine 
description, the program presented here is incomplete and represents 
at best a format for further use developing the complete solution to 
the LPCM utilizing the integral equation metheds of this thesis. 

One important result of the use of this subroutine is that the 
total computation time necessary for solution is very short, Sub- 
routine PLPBV, the MAIN program calling it, and all cf the subroutines 
used compiled under WATFOR in 12 seconds and executed completely in 
15 seconds. More elaborate programs would undoubtedly require more 


than 15 seconds computation time, but the fact that a problem as general 


fe 


as PLPBV can be solved in a time this short is very encouraging for 


more general versions of the LPCM. A flowchart of PLPBV is presented 
in Chapter L2, 

For completeness, tne MAIN vrogram used to call PLPBV as a test 
is presented at the end of this anpendix. In general, this main 
program would employ the approximations and combinations of upper and 


lower solutions for a complete distillation column modeling using 


PLPBV to solve the upper and lower versions of the LPUM, 
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“SUBROUTINE PLPBV (PRMT,P-sQyNve BC,F US UT U) 


NOTE 5 5% 6 Be XK a 2 Oe Ie eo ae 2 ae ec eae he 2 2 ae 3h ae 2 a ake ae 2h ae ae aie ae ae ae fe aie hc he ae a a kc 
THIS PROGRAM IS INCOMPLETE IN SEVERAL ASPECTS ANDO 
REQUIRES EXTENSIVE WORK TO BE USEFUL IN GENERAL 
PURPOSE COMPUTATION. 

SLBROUTINE PLPBV IS PRESENTED HERE AS A FORMAT FOR 

AND EXAMPLE QF THE STEPS NECESSARY TO APPLY THE 


INTEGRAL EQUATION TECHNIQUE OF THIS THESIS. 
NGM TISTICTSCCOOC©§OOSOCOTOSS OO OS SSCS SS PSST Se TS EL Se eer 


hI RO Ok kat ae aie ke ake a ee RI ke i a ate ae kee ake a eR a aR a a i ae ke ae i a a ak ae 2k 
SUBROUTINE PLPBV(PRMT,P.QsNoBCeFeUS,UT,U) 


PURPOSE 
TO SOLVE A LINEAR,SECCND CRDER;PARABOLIC, 
PARTIAL DIFFERENTIAL EQUATION WITH BOUNDARY 
CONDITIONS AND POLYNOMIAL COEFFICIENTS, 
D/DXX(P{XIUCX oT) 3 4D/DX( QE X UXT) P=D/DT(UCX, Td) 
WHERE P AND © ARE POLYNOMIALS IN X, IN RESPONSE 
TO STEP INPUTS»WHERE F IS STEADY STATE AT T=TL. 


USAGE 
CALL PLPBVIPRMT oP sQeNoBCoF,US, UTS) 


DESCRIPTION OF PARAMETERS 
F ~INPUT STEADY-STATE DISTRIBUTION AT TL 
BC ~INPUT BOUNDARY CONDITIONS OF THE FORM 
BC(7)UFBC(9YUX=BC(1) AT X=XuU 
BC(4)U+BC(6)UX=BC(2) AT X=XL 
BC MUST BE AT LEAST OF DIMENSION 10 
PLPBV SETS BC{(3s58,10)=0 FOR LBVP USE 
PRMT - INPUT VECTOR WHICH SPECIFIES THE PARA- 
METERS OF THE X AND T INTERVALS AND OF 
THE ACCURACY FOR SUBROUTINES USED BY 
PLPBVs MUST BE AT LEASt DIMENSION 10 
PRMT{L)J— LOWER BOUND XL OF THE X VARIABLE 
PRMT{1)- UPPER BOUND XU OF THE X VARIABLE 
PRMT{3)— SPATIAL INCREMENT OF THE X VARIABLE 
PRMT(4)— UPPER ERROR BOUND FOR X 
PRMT(5)— TERMINATION PARAMETER FOR X 
PRMT(6)-— LOWER BOUND TL OF THE T VARIABLE 
PRMT(7)- UPPER BOUND TU OF THE T VARIABLE 
PRMT(8)-— TIME INCREMENT OF THE YT VARTABLE 
PRMT(9)- ERROR WEIGHT LTE 1.0 FOR LBVP USE 
PRMT(10} ERROR WEIGHT LYE 1-0 FOR LBVP USE 
BOTH ERROR WEIGHTS ARE USUALLY 1.0 
P ~ INPUT VECTOR SPECIFYING THE 
COEFFIGTENES OF THe Erkan 
POLYNOMTAL OF DEGREE N 
Q - INPUT VECTOR SPECIFYING THE 
COEFFICIENTS OF THE SECOND 
POLYNOMIAL OF DEGREE N 
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N —- INPUT SPECIFYING THE MAXIMUM 
DEGREE OF P AND Q 
THE COEFFICIENTS MUST BE IN THE 
ORDER CF LOW TO HIGH INDEX 
CORRESPONOING TO LOW TO HIGH 
DEGREE IN EACH POLYNOMIAL 
BOTH P AND Q MUST BE AT LEAST OF 
DIMENSION N WITH THE REMAINING 
COEFFICIENTS IN EITHER P OR Q SET=0 


US —- QUTPUT STEADY-STATE SPATIAL 
DISTRIBUTION VECTOR 

UT —- QUTPUT TRANSTENT RESPONSE 

U —- OUTPUT U=US+tUT OVERALL SOLUTION TO THE 


TWO POINT BOUNDARY VALUE, STEP INITIAL 
VALUE PROBLEM. 


REMARKS 

(1) STEADY-STATE PORTION WRITTEN FOR MAX N=3, 
THIS CAN BE EASILY CHANGED BY USING 
DIFFERENT DIMENSIGN STATEMENTS ON P AND Q 

(2) TRANSIENT PORTION WRITTEN FOR P(1),Q1), 
AND O(23+4 CHANGING THIS REQUIRES ANALYTICAL 
APPLICATION OF THE INTEGRAL EQUATION 
TECHNIQUE PRESENTED IN THIS THESIS. 

(3) THIS SUBROUTINE HAS BEEN INITIALLY WRITTEN 
TO EVALUATE AT TEN POINTS IN X AND AT NINE 
POINTS IN T. THIS COULD BE CHANGED TO A 
MORE GENERAL METHOD IF DESIRED. 

(4) THERE ARE SEVERAL ERROR QUTPUTS IN PLPBV. 

(5) SUBROUTINE PLPBV HAS NOT BEEN OPTIMIZED. 

(6) SUBROUTINE PLPBV HAS BEEN WRITTEN MAINLY TO 
DEMONSTRATE THE STEPS NECESSARY TO APPLY 
THE INTEGRAL EQUATION SOLUTION OF THIS THESIS. 

(7) IT ITS EXPECTED THAT SEVERAL *BUGS* REMAIN 
IN SUBROUTINE PLPBV AS WRITTEN HERE. 

(8) AS THE AFCT OF APPENDIX AS IS NOW WRITTEN, 
A STATEMENT, CGMMON /S/P,QsN -¢IS REQUIRED 
IN THE MAIN PROGRAM. THIS COULD BE CHANGED 
BY INCREASING THE PARAMETER DIMENSIONS IN 
SUBROUTINE AFCT. ' 


SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
SUBROUTINE PLPBY REQUIRES IBM/SSP SUBROUTINES- 
LBVP. GELG. EIGENs AND LOC 
AND USER FURNISHED SUBROUTINES- 
AFCT, BDFCT. FCT, AND OUTP FOR USE IN LBVP 


2 ak ai ak ae ae age ake cake fe ake ake 2 akc aK akc ae ae ae 2k ae ae ak a ake ae ae a ae ae ake ae ake ak ake ake ae ake ae ok a a ae ake ae ak oe 2k ake ak akc oe 


DIMENSTON PRMT (1) oP01)2001),8C(1),FC1),USC1) »UT(1) .Ut1) 
DIMENSION D(4) ,AK(100)-.UZ(10),X1{100),XJ(100) 
DIMENSION AKS(55) .THE10),E1G(100),TLAM(10) swWI(100) 


INPUT AND PARAMETER SETUP SECTION 
NDIM=2 
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OO 


NV=10 
NT=9 
CN=10.0 
CT=8.0 
EPS=1.0E-05 
AA=PRMT(1) 
BB=PRMT(2) 
TL=PRMT (6) 
THH=PRMT(7) 
IF(ABS(PC1) )—EPS) 33,334.34 
33 WRITE(6.35)P(1) 
35 FORMAT(1X,*§P{1) INVALID, P(1l)= *,E£15.79//) 
RETURN 
34 SPZ=SORT(P(1)) 
SSPZ=SORT(SPZ) 


STEADY STATE SECTION. USES SUBROUTINE LBVP 
STEADY-STATE VALUES RETURNED THROUGH XPAR 
DOI NENSTON Y(2) -DERY (2) sAV14) -AUX{( 20.2) 
DIMENSTON Rf2),B8(4)5C(14),TPAR( 15) »XPAR(15) 
EXTERNAL FCT,OFCTsAFCT,OUTP 

DO 55 {Y=1,5 

XPAR(I)=PRMT(1) 

TPAR(T)=PRMT(i+5) 

55 CONTINUE 

R(1L)=BC(1) 

R{23=BC(2) 

BC(3)=0.0 

BC(5)=0.0 

BC (8)=0.0 

BC{10)=0.0 

DC 56 I[=1+44 

B(I)=BC{ f+2) 

C(1}=8C(1+6) 

56 CCNTINUE 

DERY{1L)I=PRMTI9) 

DERY{2)=PRMIT710) 

CALL LCBVPIXPAR «eBsCoReYeDERYesNDIMsIHLEF,AFCT,FCT,DFCT, 
LOUTP, AUXs AY) 


TESTS APPLIED UPON RETURN FROM LBVP 
IF (IHL F-13) 39440.41 

41 WRITE (6.42) 

42 FORMAT(1Xs'LBVP HAS IHLF=14, NO SOLUTION ',//) 
RE TURN 

39 TFCTHLF-11)44245240 

45 WRITE (6546) 

46 FCRMAT(1Xe*L BVP HAS THLF GT 16, NO SOLUTIONS,//) 
RETURN 

40 WRITE(6,43) 

43 FORMAT(1X,*LBVP HAS THLFH=H13 OR 1l2,y PARAMETER ERROR*®,//) 
RETURN 

44 WRITE(6.47) THLF 

47 FORMAT(LX,s*LBVP RETURN SATISFACTORY, IHLF= ‘%,139//) 
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Oa 


OO 


Om 


Sia Sewell akek ele ke 


oak 
32 


30 


50 
52 


D1 


TRANSTENT SECTIONs USES SUBROUTINES EIGEN AND GELG 


TRANSFORMING BOUNDARY CONDITIONS 
D€1)=B8C(7)-€BC(9)*(001)40(2) #PRMT(2)))/(2.0*P(1)) 
D(2)=BC(4)-—(BC(6) *(001)+Q(02) *PRMT(1)))/(2.0*P(1)) 
D(3)=BC(9)/SPZ 

D(4)=BC(6)/SPZ 

DTEST=0(1 )*D(4)-0(2)*D(3) 

CONE=D01)*002)-OTEST 

IFC ABST CONE J-EPS) 31-31, 30 

WRITE (6.32) CONE 

PCRMAT(1X,*°Cl IS *,E15.7e'* SETTING C1 TO 1.E-05 ",//) 
CONE=EPS 

SET UP INTERVAL FOR KERNEL EVALUATION 

DELTX={ PRMT{2)-—PRMT(1))/CN 

CC=DELTX*CN/P({1) 


PROGRAM RUNS FOR 10 DELTX AND 10 DELTZ 
FOR GREATER ACCURACY MORE INCREMENTS COULD BE USED 
DELTZ=CC/CN 


BEGIN KERNEL EVALUATIONS 
EMZ=1.5*O(2)4+06.25*Q(1)*Q(1)/Pi1) 


SEVERAL ALTERNATIVE KERNEL EVALUATION ROUTINES 
MUST BE DESIGNED AND PLACED IN THIS SECTION FGR 
CASES SUCH AS #* 

(1) M(O0)=0 

(2) MOZ)=MI0) 
ANC OTHERS. 
THESE ROUTINES MUST USE EQUATIONS Al.23 AND Al.24% 
UEP into VHESTS. 


TEST M{Q) 
IFCABS(EMZ)-EPS)56,50,51 
WRITE6e529EMZ 
FORMAT(1Xe?M(0)=0 » INVALID FOR PLPBVy MO)=%9E15.71//) 
RETURN 

DK=0(1)4D3) 

FP=4.0*P(1) 

LK=0 

i S=0 

DQ 49 JK=1,NV 

DC 48 IK=1,NV 

COMPUTE K4(Z+sS) 

LK=LK+1 

Z={IK*DEL TZ 

S=JK*DELTZ 

XS=SPZ*S+AA 

€NS=CC-S 

DMDZ= (D(4)-D{2)*Z)/CONE 
XC M¥S=SPZ*¥CMS+AA 
EMS=EMZ+{2.0¥*Q(1) *0(2)¥*XS+0(2)*OQ(2)*XS¥XS)/EFP 
XQ=0(2)*XCMS 


Pee 





MONA MANY 


54 


Zio) 


48 
49 


58 


57 
14% 


59 


60 


54 
63 


EMCMS=EMZ4+(2.0*0 (1) *XQ+XO*XQ) / FP 
KF=DMDZ* (001) *EMCMS+D( 3) *EMS—DK*EMZ)/EMZ 
IF (IK-JK)54,54,53 


AK( LK) =KF 

LS=LS+1 

AKS(LS)=KF Note: obvious error 

GQ TO 48 KF should be RKF, real. 
ZMS=Z-S 


XZMS=SPZ*ZMS+AA 

XOZ=Q(2)*XZMS 
AK(LK)=(02.0*Q(1)*XQZ4+XOZ*XOZI/(EMZ*FP) )+KE 
CONTINUE 

CONTINUE 


CCMPUTE EIGENVALUES AND EIGENVECTORS OF AKS 
MV=0 
CALL ETGENCAKS»EIGs NV oMV) 


TAKE EIGENVALUES TH(I) FROM DIAGONAL GF AKS(1) 
EPTH=1.0E-10 

[IS=] 

DO 60 f=l »NV 

THCIJ=AKSCIS) 

IS=I1S¢( 141) 

IFC THOT )-EPTHI57,57-58 
TLAM(T)=1.0/( THC I)*DELTZ) 

GO TO 60 

WRITE(6, 74) THC I),.1 

FORMATC1Xe* THETACI) INVALID, TH(I)="',£15.7,° l=*, 13) 
WRITE(6,59) 

FORMATCLX,* SETTING THI IT)=EPTH' ,/) 

TLAM(TJ=EPTH 

CONTINUE 


THIS PROGRAM HAS BEEN SET UP SO THAT DELTX AND DELTZ 
CORRESPOND. THIS MAKES WI{ZTJHWICXT). IF THIS WERE 
NOT DESIRED, THeN A SECTION UTILIZING EQUATION A3.10 
WCULD HAVE TO BE EMPLOYED. 


CCMPUTE WI@X) FROM WI(Z)ie THETA DI=THOT) WI(Z)=EIG(Z) 
IJ=0 

DO 63 T=lsNV 

DO 64 J=1s»NV 

[J=fJ+1 

WITT JISETGCIJ) 

CONTINUE 

CONTINUE 


FRANSFORM WI(X) TO XIX) 

X=DELTX+AA 

DC 6&1 J=lsNV 

XMA=X—-AA 

XP A=X+AA 

XFUN=( XMA*Q(1)9/02.08PC1) +4 XMA*¥XPA*Q(2))/FP 


EXX=EXP(—XFUN) ¥SSPZ 
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mon 


CVO? CHES 


ap I a ee 1 «a 


62 
61 


65 


67 


36 
68 


38 
69 


37 


70 


X=X+DELTX 

DO 62 IT=1.NV 
JI=J+Nv*(I-1) 
XICSLTISEXX¥WIC ST) 
CONTINUE 

CONTINUE 


BEGIN EIGENFUNCTION EXPANSION TO MEET INITIAL CONDITION. 
DELTT=({ THH-TL)/CT 
T=TL 


FORM STEADY-STATE VECTOR 
DO 65 [T=1,NV 
UZCTI=FCT)-US(T) 
CONTINUE 


FORM EXPANSION MATRIXsXJ(X)_, AT T=TL 
IJ=0 

DO 66 I=1.NV 

EXT=EXP(-TLAM( 1) ¥*T) 

DC 67 J=leNV 

TJ=fJ+l 

XJCISI=X ICT S* EXT 

CONTINUE 

CONTINUE 


SOLVE SIMULTANEOGUS EQUATIONS FOR EIGENFUNCTION 
CCNSTANTS USING SUBROUTINE GELG 

NG=1 

CALL GELG{UZ»sXJoNVeNGeEPS,IER) 


TEST ITER UPON RETURN FROM GELG 

IFC IER 136537238 

WRITE(6,68) 

FORMAT(1X,"*NO RESULTsPIVOT ELEMENT=0O IN GELG *,/) 
RETURN 

WRITE 16,69} 

FORMAT (LXse* POSSIBLE LOSS OF SIGNIFICANCE IN GELG*,/) 


TRIS SECTIGN SETS UP UC(X.T) IN A FORM ACCEPTIBLE 
FOR PLOTTING BY SUBROUTINE PLOT QF APPENDIX. A4. 
TRE FIRST COLUMN OF UWI) TS THE INDEPENDENT 
VARIABLE Xs THE SECOND COLUMN IS THE STEADY-STATE-s 
AND THE REMAINING 8 COLUMNS ARE RESPONSES. 

X=AA 

DC 70 [=1,NV 

K=X+DELTX 

UTI )=X 

UC I+1lO3=F (1) 

CONTINUE 


CALCULATE TRANSIENT SOLUTION MATRIX UT{XeT) 


SET TIME-SPACE GRID 
I JK=20 
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OAD 


Lae 


12 
v1 


T=TL 

NTC=NT-1 

OC 71 I[=1eNTC 
T=T+DELTT 
X=AA 

DO 72 J=leNV 
X=X+DELTX 


TIME-SPACE GRID NOW SET 


SUM SERIES FOR EACH VALUE 
IJK=[SK+tl 

UT (IJK )=0.0 

DO 73 K=1+NV 

I X=J+NV*CUK-1) 
EXT=EXP(-TLAM(K) *T) 

UTC I SKI=UZ(K)* XI CTX) FEXT+HUT (CI UK) 
CONTINUE 
UCTJKI=US(JIFUTCI SK) 
CONTINUE 

CONTINUE 

RETURN 

END 
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$ JOB 
C 
€ 
C 
C 
100 
101 
102 
C 
‘ 
10 
c 
¢ 


HAYES »sKP=29.TIME=2,PAGES=50 


MAIN PROGRAM TQ USE PLPBV AS A TEST 
DIMENSTON BC{(10).U(100),US(10),UT(100) ,F( 10) 
DIMENSTON P(3) .0(3) eAl2+4).PRMT (20) 

CCMMON /S/P.0. M 


PLPBV INPUT SECTION 
M=2 
READ{5.100)(P(I) O11) »IT=1)M) 
FORMAT(2F10.5) 
READ(5+101) ((A0CT.J5) oe J=104) 51T=1,2) 
FORMAT (4F10.5) 
READ(5.102) (PRMT(I).f=1.10) 
FORMAT(2E15.7) 
BC(7)=AC1,51) 
BC (9) =A(1,.2) 
BC(1)=A(1.4) 
BC(4)=A(2,13 
BC(6)=Al2,2) 
BC{(2)=A(2.4) 


SET UP LINEAR STEADY-STATE AS A TEST 
X=0.0 

00 10 KSS=1,10 

F(KSS)=0.0 

X=X+0.1 

US(KSS)=X 

CONTINUE 


TEST STEADY-STATE INPUT FOR PLPBV FROM UC1L,APPENDIX AG 
PRMT(11)=0.06 | 
PRMT{12)=0.126 

PRMT(13)=0.200 

PRMT@14)=0.27 

PRMT{15)=0. 37 

PRMT{16)=0.47 

PRMT(17)=0.59 

PRM7{118)=0.71 

PRMT(19)=0. 84 

PRMT{20)=1.00 

CALL PLPBV( PRMT oP 5QeM,BC,FeUSeUT,U) 

CALL EXIT 

END 
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SECTION 6 


BIBLIOGRAPHY AND REFERENCING (B) 


Bl BIBLIOGRAPHY 
B2 BIBLIOGRAPHY REFERENCES LISTED BY AREA OF APPLICATION 


B3 BIOGRAPHICAL NOTE 


AN UNDERSTATEMENT: 
“ THE LITERATURE ON DISTILLATION IS VOLUMINOUS - " 


R. Je HENGSTEBECK (H-18) 


ONE OF THE MAJOR EFFORTS OF THIS THESIS TURNED OUT TO BE THE 
COMPILATION OF THE BIBLIOGRAPHY PRESZNTED IN THIS SECTION, THE 
REFERENCES LISTED APPLY TO FOUR MAJOR AREAS PERTINENT TO THIS 
THESIS. 

1. GENERAL THEORY OF DISTILLATION 

2. DISTILLATION COLUMN DYNAMICS 

3. DISTILLATION COLUMN CONTROL 

4, MATHEMATICS AND COMPUTATION 
THE AUTHOR HOPES THAT UTILIZATION OF CHAPTER B2, BIBLIOGRAPHY 
REFERENCES LISTED BY AREA OF APPLICATION, WILL RESULT IN CONSIDER~ 
ABLE SAVINGS OF SEARCHING TIME AND &FRPORT FOR ANYONE INTERESTED IN 


THESE AREAS, 
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CHAPTER Bl 


BIBLIOGRAPHY 

This chapter presents a bibliography of 352 references pertinent 
to the general area of distillation and to the specific areas of this 
thesis. The references consist only of those in English, although 
extensive literature on distillation has been published in the foreign 
journals, especially in Russian and German. For the most part, the 
references were taken from the following journals for the years from 
about 1955 to 1969. 

1. Industrial and Engineering Chemistry 

2. American Institute of Chemical Engineers 

3. Chemical Engineering Progress (and Symposium Series) 

4, Transactions of the Institution of Chemical Engineers 

5. Chemical Engineering Science 

6. British Chemical Engineering 

7- Canadian Journal of Chemical Engineering 

The entries in this bibliography are in the order of the first 
letter of the author's last name, for the reader's mnemonic convenience, 
but no attempt has been made to subalphabetize within each group for 
this author's convenience, 

It is this author's intention that Chapter B2 - Bibliography 
References Listed by Area of Application be used in conjunction with 
this chapter for any given subject or area of research. In each ref- 
erence in this bibliography an attempt has been made, especially with 


journal articles, to present as complete a description as possible of 
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the reference, even at the cost of complete overspecification and 
possible redundancy of information. 

The criterion of availability of each reference in the M.I.T. 
libraries (except for theses) placed a significant constraint on the 
number of references which have been listed here. The author chose 


this as his stopping point in the compilation of this bibliography. 
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CHAPTER B2 


BIBLIOGRAPHY REFERENCES LISTED BY AREA OF APPLICATION 
The purpose of this chapter is to classify the references in 
Chapter Bl, the Bibliography, according to their principal areas of 
application. The general format of this classification is presented 
in Table B2.1. Many of the references apply to more than one area and 
are so listed for the user's convenience, Rather than merely list the 
references under each area in alphabetical order, the author has chosen 
to present them column by column in the order of decreasing utility to 
the study of the subject or of decreasing clarity. In other words, the 
first listed in each area seem to the author to be the most important 
for anyone researching that area. Recognizing that such a listing is, 
indeed, very subjective, the author apologizes to anyone who may find 
them “out of order" with respect to his particular slant on the subject. 
B2.1 General Theory of Distillation 

Textbooks 

Extensive Bibliographies and Literature Surveys 

References of Historical Interest 

General Distillation 
Dynamic or Transient Analyses 

Steady-State Analysis and McCabe-Thiele Diagrams 

Structural Design 

Economics and Operations Analysis 

Thermodynamics 

Hydrodynamics 


Chemistry 


Philosophy Table B2.1 (Contd. ) 








B2.2 Distillation Column Dynamics 


Textbooks 
Theses 
Reviews, Bibliographies, and Literature Surveys 
Dynamic Models or Solutions 
Discrete Plate Equations 
Frequency Analysis or Laplace Transform Solution 
Transient or Time Analysis 
Numerical Solution 
Digital Computer 
Hybrid Computer 
Analog Computer 
Analytical Analysis 
Continuous Spatial Equations 
Frequency Analyses or Laplace Transform Solution 
Transient or Time Analysis 
Experimental Transient Behavior 
Frequency Response 
Time Response 
vere Distillation 
B2.3 Distillation Column Control 
Textbooks 
Theses 
Extensive Bibliographies and Literature Surveys 
Conventional Control Systems 
Digital Control 
Hybrid Control 


Table B2.1 (Contd. ) 
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Analog Control 
Instrumentation 
Control Systems Using Dynamic Models 
Digital Computer Control 
Hybrid Computer Control 
Analog Computer Control 
Optimal Control 
Distributed or Modal Control 
B2.4 Mathematics and Computation 
Ordinary Differential Equation Theory 
Partial Differential Equation Theory 
Integral Equation Theory 
Mathematical Transformations 
Matrix Mathematics 
Numerical Solution Techniques 
Boundary Value Problems 
Higen - Values, Vectors, and Functions 
special Functions 


Computation and Computer Programming 


Table B2.1 = FORMAT FOR CLASSIFICATION BY AREA OF APPLICATION 
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B2.1 GENERAL THEORY OF DISTILLATION 
Textbooks 
(B-1)  (A-19) (H-18) (R-22) (T-1) 
(G-3)  (B-15) (L-26) (R-17) (B-12) 
(H-10) (M-14) (C-5) (S-8)  (H-5) 
(R-23) (0-2) (C-1) (8-9)  (C-7) 
(R-21) (H-2) (P-3) (S-3)  (M-2) 
(V-2) (H-9) (P-10) (S-1) (V¥-3) 
Extensive Bibliographies and Literature Surveys 
(W-16) (G-9)  (W-13) (W-18) (F-6) 
(B-17) (R-12) (W-15) (W-26) (G-14) 
(B-18) (P11) (F-7) (Z-3)  (W-23) 
(B-19) (R-8) (R-25) (W-8) 
(W-10) (H-7)  (R-18) (W-17) 
(W-11) (W-9)  (H-8)  (G-10) 
(F-8) (W-12) (R-17) (wW-19) 
References of Historical Interest 
General Distillation 
(L-19) (L-25) (V-2) (2-3) (P12) (G-14) (A-17) 
(M-10) (M-9)  (R~29) (C-15) (T-6) (I-1) (U-1) 
Dynamic or Transient Analysis 
(M-8)  (B-29) (L-2) (B-7) (R-2) (R-12) (J-1) (P-5) 
Steady State Analysis and McCabe - Thiele Diagrams 
(B-1)  (H-22) (C-1) (R-21) (8-23) (L-25) (S-7)  (R-18) 
(M-10) (B-10) (F=3)  (H-30) (S-12) (R-30) (S-29) (A-19) 
(E-4)  (M-18) (F-4) (3-7) 9 (C-14) (T-1) — (W-30) = (P-12) 


(H-10) (P-8) (V-2) (L-12) (D-2) (H-2) (F-9) (T-7) 
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B2ec 


Steady State Analysis and McCabe - Thiele Di 


(H-18) (V-2) (F-15) (A-10) 
(H-26) (S-15) (G-3)  (A-17) 
(C-9)  (S-20) (H-5)  (B-24) 
(M-8) (2-2) (H-2) (B-26) 
29 eB=26) 


Structural Design 
(L-13) (R-22) (F-10) (J-3) 


Bconomics and Operations Analysis 


(L-15) (M-13) (B-15) (S-8) 

Thermodynamics 
(H-9)  (R-23) 

Hydrodynamics or Fluid Mechanics 
(F-4)  (V-2)  (R-23) (S-5) 
(L-26) (B-12) (H-14) (G-121) 

Chemistry 
(B~16) (B-1) (H-17) (T-6) 
(H-8)  (H-10) (F-9)  (G-8) 

Philosophy 


(A-13) (B-28) 


DISTILLATION COLUMN DYNAMICS 


Textbooks 
(M-1)  (H-7) (G-3) = (F-16) 
Theses 
(F-14) (W-21) (D-14) (C-10) 
(M-15) (M-12) (M-17) (Q-1) 
(S-28) (A-5)  (R-6) = (R-27) 
(S-17) (C-4) (A-16) (L-17) 
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(0-2) 
(T-6) 
(P=13) 
(R-5) 


(M-11) 


UTS), 


(P-9) 
(F-5) 


(S-5) 
(H-20) 


(B-32) 
(0-3) 
(S-22) 


($-18) 


ams (Contd. 
(S-15) (F-1) 
(S-5) 
(S25) 
(S-24) 


(H-21) 


(L-12) (D-11) 


(F-10) (C-7) 


(B-4) - (S-19) 


(M-24) (T-7) 
(R-23) 
(T-6) 


(W-22) (S-27) 
(G-5)  (W-20) 
(M-19) 


(M-24 ) 


(W-24) 


(I-1) GEG 
(L-24) (S-1) 





Reviews, Bibliographies, and Literature Surveys 


(A-3)  (H-7) = (R25) 
(R-12) (B-17) (R-8) 
(W-16) (B18) (W-10) 


Dynamic Models or Solutions 


Discrete Plate Equations 


(W-15) (Z-3) (S-33) (L-3) 
(W-12) (G-9) (W-13) (T-2) 
(W-11) (B-19) (W-24) (W-23) 


Frequency Analysis or Laplace Transform Solution 
(R<25) (W-15) (Z-3) (S-33) (L-3) 
(R-12) (B17) (R-8) 


(A-3)  (H-7) 


(W-16) (B-18) 


(W-12) (G-9)  (W-13) (T-2) 


(W-10) (W-11) (B-19) (W-24 (W-23) 


Transient or Time Analysis 


Numerical Solution 


Digital Computer 


(D-4) 
(B-32) 
(R-2) 
(D-8) 
(H-3) 
(L-5) 
(R=-20) 


(D-15) 


(R-16) 
(R-9) 
(A-4) 
(M-3) 
(R-10) 
(R-11) 
2) 


(W-1) 


Hybrid Computer 


(72H) 


(F-13) 


Analog Computer 


ERE S) 


(P-1) 
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(P-2) (T-9) (S16) 
(R-17) (W-25) (R-19) 
(De9) (¥-1) = (L~27) 
(L-18) (P-5) (S-20) 
(D-6) (G-4) — (W-4) 

(A-9)  (R-5) = (R18) 


(S-6)  (D-11) 
(S-33) (R-1) 
(R-6) (F-12) 


(R-3)  (G-12) (L-5) 





2.3 


Analytical Analysis 
(D-3)  (M-1)  (S-31) (B-27) (H-20) 
(G-2) (R=7) (W-2) (C-2) = (H-21) 
(A-3)  (R-8) = (F-16) (2-3) = (Be?) 
(M-8)  (R-32) (G-3) (F+5)  (B-12) 
Continuous Spatial Equations 
Frequency Analysis or Laplace Transform Solution 
(5-6) (Me27) (Het) (s+?) (We5) (+2) 
Transient or Time Analysis 
(J-1) (R-8) 9 (Mel) = (L-26) (T-5) — (J-5) (W214) 
(H+)  (D=3) (M-28) (P=7) (G-3) (L-20) (B-25) 
(K-2)  (S-13) (K=-3) (R-14) (H-25) (B-20) (C-8) 
(O-1) (D-14) (C-2) (Be29) (J-2) (B-24) 
Experimental Transient Behavior 
Frequency Response 
(H-15) (H-23) (A-6) (H-24) (W-28) 
TE Response 
(H-3)  (L-5) (Be5) (Be3) 9 (R-15)  (S-4) 9 (D-12) — (R-30) 
Cyclic Distillation 
(M-19) (A-14) (M-20) (B-17) (B-18) (S-4) 
DISTILLATION COLUMN CONTROL 
Textbooks 


(A-12) (G-3) (Be13) (K=5) (C-1) (C#15) (A-211) (L-1) 


Theses 


(B-2) (Be32) (M-5) (G-1) (S-18) (M-25) 
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Extensive Bibliographies and Literature Surveys 


(A-3)  (B-17) (R-25) (R-12) (W-10) (W-20) (W-12) (W-21) 


(W-16) (B-18) (R-8) (S-33) (W-11) (B-19) (W-15) (W-23) 


Conventional Control Systems 


Digital Control 
(S-32) (H-29) (B-6) (L-15) (L-24) (7-8) 

Hybrid Control 
(F-13) (F-1) (F-12) 

Analog Control and Instrumentation 
(H-16) (L-16) (B-10) (C=15) (G-12) (S-30) (M-29) (S-9) 
(R-31) (B-13) (B-31) (P-6) (K-4) (C-12) (P-14) (W-23) 
(G-3)  (B-22) (C-1) (S-1) (L-4) (T-11) (R-13) (H-1) 


Control Systems Using Dynamic Models 
Digital Computer Control 

(A-12) (D-4)  (R-20) (L-18) (D-6) (S=-33) (C-11) (¥-18) 
Hybrid Computer Control 

(D-13) (A-8) 
Analog Computer Control 

(J-6) (L-14) (L-6) (W-27) (H-15) (K-2) (L-3) (R-3) 
Optimal Control 

(L-9)  (B-32) (D-13) (A-8)  (S-33) (J-5) 

(B-11) (A-11) (K-5)  (w-29) 
Distributed or Modal Control 


(D-9) (J-6) (M-5) (S-31) (C-4) (S-2) (F-16) (G-13) 
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B2.4 MATHEMATICS AND COMPUTATION 


Ordinary Dirrerential Equation Theory 
(H=-11) (S-21) (B-14) (I-2) (H-12) 


Partial Differential Equation Theory 


(G-7) (B-21) (0-1) (F-16) 


(F-11) (A-15) (F-12) (G-6) 
Integral kquation Theory 
(T-3)  (L-8) (D-5) 
(P-4)  (H-13) (M6) 


Mathematical Transformations 


(S-11) 
(S-14) 
(T-4) (2-1) (S-31) (M-16) 
Matrix Mathematics 

(A-4)  (A-12) (H-13) (J-7) 
Numerical Solution Techniques 

(B-30) (F-11) (M-22) (P-5) 
(D-7)  (A-10) 
(S-13) 


(C-8) 


(A-15) (K-1) 


(R-24) (L-10) (R=l4) 
(D-8) 


(H-28) 


(L-1) (P-8) 
(R-29) 


(P-12) 


(H-19) 
(M-4) 


(L-21) 
(F-2)  (B-32) 
Boundary Value Problems 
(L-11) (B-8) (S-10) (H-12) 
(H-25) 


(@e=22) 


(H=-11) (B-14) (V-1) 


(K-1)  (B-32) (F-2) 
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(W-7) 
(H-28 ) 


(V-4) 


(S-2) 


(R-14) 


(M-26) 


(H=25) 
(J-7) 
(L-11) 
(L=22) 
(L-7) 
(A-17) 


(L-7) 
(B-9) 
(C-6) 


(S-10) 


(T=10) 


(G-7) 


(C-8) 


(D-9) 
(D-10) 
(M-2) 
cy) 
(R-17) 
(T-10) 


(M-22) 


23) 


(H-6) 


(W-6) 
(W-7) 


(W-26) 
(M-26) 
(J-5) 
(D-15) 
(S-12) 











Eigen - Values, Vectors, and Functions 
(B-30) (H-11) (B-14) (m-2)  (S-32) 
(A-12) (H-13) (C-13) (8-13) (v-1) 
(G-13) (1-2) (D-3) (S-10) (w-26) 

Special Functions 
(A-18) (L-23) (M-16) (R-28) (H-27) 

Computation and Computer Programming 
(M-4)  (I-3)  (C-13) (M-18) (B15) 
(O-4)  (P-8)  (H-29) (S-32) (M-25) 
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(C-6)  (J-7) 
(B-5) 
(H-12) 
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CHAPTER B3 
BIOGRAPHICAL NOTE 

Lt. Michael N, Hayes was born in Abilene, Texas, on September 2, 
1941, Four years later the family moved to Ohio, where he attended 
primary schools in Wakeman and in Massillon, He was graduated from 
Washington High School, Massillon, Ohio, in June 1959. 

Mr. Hayes enlisted as a Seaman Recruit in the United States Navy 
in September, 1959. He was graduated from the U.S. Naval Guided 
Missiles School, Dam Neck, Virginia, in July, 1960, and was assigned 
to the U.S.S,. Kitty Hawk (CVA-63) until June, 1961. He was graduated 
from the U.S. Naval Preparatory School, Bainbridge, Maryland, in 
August, 1961. 

In September, 1961, Mr. Hayes began 4 years of undergraduate 
education at North Carolina State University, Raleigh, N.C., under 
the auspices of the Naval Enlisted Scientific Education Program 
(NESEP), He was graduated in June, 1965, with the degrees of 
Bachelor of Science in both Electrical Engineering and Applied 
Mathematics and with the degree of Master of Electrical Engineering. 

After attending the U.S. Naval Officer Candidate School at 
Newport, R.I., Mr. Hayes was commissioned an Ensign in the U.S. Navy 
in October, 1965, and assigned for duty at the Boston Naval Shipyard, 
Boston, Massachusetts. While stationed at the Shipyard, he worked 
as a Ship Superintendent, supervising naval ship construction and 
repair, and attended evening classes at Boston University in business 


administration and at Northeastern University in pure mathematics. 
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In June, 1967, Mr. Hayes entered M.I.T. as a U.S. Naval Post- 
graduate student in Curriculum 13A with specialized studies in 
electrical engineering. He is a member of Tau Beta Pi, Eta Kappa Nu, 
Phi Eta Sigma, and is an Eagle Scout, 

Lt. Hayes is married to the former Ann Marie Graney of Stoughton, 
Massachusetts. They have a one-year-old son, Stephen, and are 


expecting their second child in December, 1969, 
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